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ABSTRACT 

Inverse-Compton cascades initiated by energetic gamma rays (E > lOOGeV) enhance the GeV emission from 
bright, extragalactic TeV sources. The absence of this emission from bright TeV blazars has been used to constrain 
the intergalactic magnetic field (IGMF), and the stringent limits placed upon the unresolved extragalactic gamma- 
ray background (EGRB) by Fermi has been used to argue against a large number of such objects at high redshifts. 
However, these are predicated upon the assumption that inverse-Compton scattering is the primary energy-loss 
mechanism for the ultra-relativistic pairs produced by the annihilation of the energetic gamma rays on extragalactic 
background light photons. Here we show that for sufficiently bright TeV sources (isotropic-equivalent luminosities 
> 10^^ ergs"') plasma beam instabilities, specifically the "oblique" instability, present a plausible mechanism by 
which the energy of these pairs can be dissipated locally, heating the intergalactic medium. Since these instabilities 
typically grow on timescales short in comparison to the inverse-Compton cooling rate, they necessarily suppress 
the inverse-Compton cascades. As a consequence, this places a severe constraint upon efforts to limit the IGMF 
from the lack of a discernible GeV bump in TeV sources. Similarly, it considerably weakens the Fermi limits 
upon the evolution of blazar populations. Specifically, we construct a TeV-blazar luminosity function from those 
objects presently observed and find that it is very well described by the quasar luminosity function at z ^ 0.1, 
shifted to lower luminosities and number densities, suggesting that both classes of sources are regulated by similar 
processes. Extending this relationship to higher redshifts, we show that the magnitude and shape of the EGRB 
above ^ 10 GeV is naturally reproduced with this particular example of a rapidly evolving TeV-blazar luminosity 
function. 

Subject headings: BL Lacertae objects: general - gamma rays: general - instabilities - magnetic fields - plasmas 
- radiative mechanisms: non-thermal 



1. INTRODUCTION 

Imaging atmospheric Cerenkov telescopes, such as H.E.S.S., 
VERITAS, and MAGIC,' have opened the very-high energy 
gamma-ray (VHEGR, E > 100 GeV) sky, finding a Universe 
populated by a variety of energetic, VHEGR sources. While 
the majority of observed VHEGR sources are Galactic in origin 
(e.g., supernova remnants, etc.), the extragalactic contribution 
is dominated by a subset of blazars. There are presently 46 ex- 
tragalactic TeV sources known^, of which 28 have well defined 
spectral energy distributions (SEDs), and are collected in Table 
1. Of these 28 well-studied objects, 24 are blazars, implying 
that blazars make up an overwhelming majority of the bright 
VHEGR sources. 

All of the extragalactic VHEGR emitters are relatively 
nearby, with z < 0.5 generally, and z^O.l typical. This is a re- 
sult of the large opacity of the Universe to TeV photons, which 
annihilate upon soft photons in the extragalactic background 
light (EBL), producing pairs (see, e.g., Gould & Schreder 
1967; Salamon & Stecker 1998; Neronov & Semikoz 2009). 
Typical mean free paths of VHEGRs range from 30Mpc to 
I Gpc depending upon gamma-ray energy and source redshift, 

' High Energy Stereoscopic System, Major Atmospheric Gamma Imaging 
Cerenkov Telescope, Very Energetic Radiation Imaging Telescope Array Sys- 
tem. 

^ See http://www.mppmu.mpg.de/~rwagner/sources/ for an up-to-date list. 



and thus the absence of high-redshift VHEGR sources is not 
unexpected. 

The pairs produced by VHEGR annihilation are necessarily 
ultrarelativistic, with typical Lorentz factors of 10^-10^. The 
standard assumption is that these pairs lose energy almost ex- 
clusively through inverse-Compton scattering the cosmic mi- 
crowave background (CMB) and EBL. When the up-scattered 
gamma-ray is itself a VHEGR the process repeats, creating a 
second generation of pairs and up-scattering additional pho- 
tons. The result is an inverse-Compton cascade (ICC) deposit- 
ing the energy of the original VHEGR in gamma-rays with en- 
ergies < 100 GeV. This places the ICC gamma rays in the LAT 
bands of Fermi, and thus Fermi has played a central role in 
constraining the VHEGR emission of high-redshift blazars. 

Based upon Fermi observations of z ~ O.I TeV sources, a 
number of authors have now published estimated lower bounds 
upon the intergalactic magnetic field (IGMF; see, e.g., Neronov 
& Vovk 2010; Tavecchio et al. 2010, 201 1; Dermer et al. 201 1; 
Taylor et al. 2011; Takahashi et al. 2012; Dolag et al. 2011). 
Typical numbers range from 10-'''G to 10"'^ G, with the latter 
values being of astrophysical interest in the context of the for- 
mation of galactic fields^. These limits on the IGMF arise from 
the lack of the GeV bump associated with the ICC of the blazar 

^ After contraction and a handful of windings, nG field strengths can be 
produced from an IGMF of ~ 10"'^ G. 
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Table 1 

List of TeV Sources with Measured Spectral Properties in Decreasing lOOGeV-lOTeV Flux Order 



Name 




uc 


JO 




a ^ 


F ^ 




a s 


At 


Class ' 


Rc f c rc nc G 


Mkn421 


0.030 


129 


68 


1 


3.32 


1.7 X 10' 


45.6 


3.15 


2.8 X 102 


H 


Chandra et al. (2010) 


lES 1959+650 


0.047 


201 


78 


1 


3.18 


1.6 X lO' 


45.9 


2.90 


85 


H 


Aharonian et al. (2003) 


lES 2344+514 


0.044 


190 


120 


0.5 


2.95 


2.3 X 10^ 


45.0 


2.82 


5.0 X 10^ 


H 


Albert et al. (2007c) 


Mkn 501 J 


0.034 


150 


8.7 


1 


2.58 


85 


44.4 


2.39 


1.6 X 10^ 


H 


Huang et al. (2009) 


3C 279 


0.536 


2000 


520 


0.2 


4.11 


68 


46.9 


2.53 


2.0 


Q 


MAGIC Collaboration et al. (2008a) 


PKS 2155-304 


0.116 


490 


1.81 


1 


3.53 


64 


45.4 


2.75 


3.3 X 10^ 


H 


HESS Collaboration et al. (2010) 


PG 1553+113 


>0.09 


> 380 


46.8 


0.3 


4.46 


41 


>44.9 


<4.29 


< 5.7 X 10' 


H 


Aharonian et al. (2008b) 


W Comae 


0.102 


430 


20 


0.4 


3.68 


31 


44.9 


3.41 


1.3 X lO'' 


I 


Acciari et al. (2009b) 


3C 66A 


0.444 


1700 


40 


0.3 


4.1 


28 


46.3 


2A3 


13 


I 


Acciari et al. (2009c) 


lES 1011+496 


0.212 


870 


200 


0.2 


4 


26 


45.5 


3.66 


2.5 X 10' 


H 


Albert et al. (2007b) 


lES 1218+304 J 


0.182 


750 


11.5 


0.5 


3.07 


24 


45.4 


2.37 


1.0 X 10^ 


H 


Acciari et al. (2010a) 


Mkn 180 


0.045 


190 


45 


0.3 


3.25 


20 


44.0 


3.17 


8.2 X 10^ 


H 


Albert et al. (2006) 


IH 1426+428 


0.129 


540 


2 


1 


2.6 


20 


45.0 


1.71 


2.1 X 10^ 


H 


Aharonian et al. (2002) 


RGB J07 10+591 J 


0.125 


520 


1.36 


1 


2.69 


15 


44.8 


1.83 


3.5 X 10' 


H 


Acciari et al. (2010c) 


lES 0806+524 


0.138 


580 


6.8 


0.4 


3.6 


10 


44.7 


3.21 


1.4 X 10^ 


H 


Acciari et al. (2009a) 


RGB J0152+017J 


0.080 


340 


0.57 


1 


2.95 


8.5 


44.1 


2.45 


3.0 X 10^ 


H 


Aharonian et al. (2008a) 


lES 1101-232 J 


0.186 


770 


0.56 


1 


2.94 


8.2 


44.9 


1.50 


2.3 X 10^ 


H 


Aharonian et al. (2007a) 


lES 0347-121 J 


0.185 


770 


0.45 


1 


3.1 


8.2 


44.9 


1.67 


2.9 X 10^ 


H 


Aharonian et al. (2007b) 


IC 310 


0.019 


83 


1.1 


1 


2.0 


8.1 


42.8 


1.90 


4.8 X lO'* 


H 


Aleksic et al. (2010) 


PKS 2005-489 


0.071 


300 


0.1 


1 


4.0 


8.0 


44.0 


3.56 


2.3 X lO'* 


H 


Aharonian et al. (2005) 


MAGIC J0223+430 






17.4 


0.3 


3.1 


7.6 




< 3.1 




R 


Aliu et al. (2009) 


lES 0229+200 J 


0.140 


590 


0.7 


1 


2.5 


6.4 


44.5 


1.51 


4.7 X 10^ 


H 


Aharonian et al. (2007c) 


PKS 1424+240 


<0.66 


< 2400 


51 


0.2 


3.8 


6.3 


<46.1 


> 1.42 


4.0 


I 


Acciari et al. (2010b) 


M87 


0.0044 


19 


0.74 


1 


2.31 


5.9 


41.4 


2.29 


1.5 X 10* 


R 


Acciari et al. (2008) 


BL Lacertae 


0.069 


290 


0.3 


1 


3.09 


5.4 


43.8 


2.67 


8.4 X 10^ 


L 


Albert et al. (2007a) 


H 2356-309 


0.165 


690 


0.3 


1 


3.09 


5.4 


44.6 


1.86 


6.5 X 10^ 


H 


Aharonian et al. (2006) 


PKS 0548-322 J 


0.069 


290 


0.3 


1 


2.86 


4.0 


43.7 


2.44 


8.4 X 10^ 


H 


Aharonian et al. (2010) 


Centaurus A 


0.0028 


12 


0.245 


1 


2.73 


2.8 


40.7 


2.72 


1.2 X 10' 


R 


Raue et al. (2010) 



and radio galaxies of Faranoff-Riley Type I (FR I), respectively. 



Comoving distance in units of Mpc 
^ Normalization of the observed photon spectrum that we assume to be of the form dN/dE = /q{£/£q)~'^, in units of 10~ 

Energy at which we normalize the spectrum, in units of TeV 

Observed spectral index at 
^ Integrated energy flux between 100 GeV and 10 TeV, in units of 10"^'^ ergcm~^s~^ 
^ Inferred isotropic integrated luminosity between lOOGeVand lOTeV, in units of ergs~^ 
^ Inferred intrinsic spectral index at Eq 

^ Time delay after which plasma beam instabilities dominate invcrsc-Compton cooling, in units of yr 
' H, I, L, Q, and R correspond to high-energy, intermediate- energy, low-energy pe;iked BL Lacs, flat spectrum radio quasars. 
Used to place limits upon the IGMF 

TeV emission, presumably due to the resulting pairs being de- 
flected significantly under the action of the IGMF itself. The 
wide range in the estimates upon the minimum IGMF is due 
primarily to different assumptions about the TeV blazar duty 
cycle. 

Fermi has also provided the most precise estimate of the 
unresolved extragalactic gamma-ray background (EGRB) for 
energies between 200 Me V and 100 GeV. Since ICCs repro- 
cess the VHEGR emission of distant sources into this band, 
this has been used to constrain the evolution of the luminos- 
ity density of VHEGR sources (see, e.g., Narumoto & Totani 
2006; Kneiske & Mannheim 2008; Inoue & Totani 2009; Ven- 
ters 2010). Generally, it has been found that these cannot 
have exhibited the dramatic rise in numbers by z ~ 1-2 seen 
in the quasar distribution. That is, the comoving number of 
blazars must have remained essentially fixed, at odds with both 
the physical picture underlying these systems and with the ob- 
served evolution of similarly accreting systems, i.e., quasars. 

Both of these conclusions depend critically upon ICCs dom- 
inating the evolution of the ultra-relativistic pairs. However, 
as we will show, the pairs constitute a cold, highly collimated 
plasma beam moving through a dense, stationary background, 
both of which are susceptible to collective plasma phenomena. 
Such beams are notoriously unstable; for instance, equal den- 



sity beams typically lose a significant fraction of their energy 
after propagating distances measured in plasma skin depths of 
the background plasma. If the VHEGR-generated pairs suf- 
fer a similar fate while propagating through the intergalactic 
medium (IGM), the cooling of the pairs would be dominated 
by plasma instabilities, thereby quenching the ICCs. 

Here we present a plausible alternative mechanism by which 
the energy in the ultra-relativistic pairs can be extracted. While 
a variety of potential plasma beam instabilities exist, we find 
that the most relevant for the VHEGR-produced pair beams is 
the "oblique" instabihty (Bret et al. 2004, 2005a; Bret 2009; 
Bret et al. 2010b; Lemoine & Pelletier 2010). This is a more 
virulent cousin of the commonly discussed Weibel and two 
stream instabilities. 

In Section 2 we discuss the formation and properties of the 
ultra-relativistic pair beam, including limits upon its temper- 
ature and density. Section 3 presents growth rates for a va- 
riety of plasma instabilities, including the oblique instability. 
Particular attention is paid to if and when plasma instabilities 
dominate inverse Compton as a means to dissipate the kinetic 
energy of the pairs. The resulting implications for studies of 
the IGMF are described in Section 4, including how such ef- 
forts might mitigate the systematic uncertainties arising from 
plasma cooling. The evolution of the luminosity function of 
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VHEGR-emitting blazars is discussed in Section 5, in which 
we construct a blazar luminosity function based upon that of 
quasars which is consistent with current TeV source popula- 
tions and the Fermi estimates of the EGRB. Finally, conclu- 
sions are contained in Section 6. 

This is the first in a series of three papers that discuss the po- 
tential cosmological impact of the TeV emission from blazars. 
Here we provide a plausible mechanism for the local dissipa- 
tion of the VHEGR luminosity of bright gamma-ray sources. 
In addition to the particular consequences this has for stud- 
ies of high-energy gamma-ray phenomenology, it also provides 
an novel heating process within the IGM. Chang et al. (2012, 
hereafter Paper II) estimates the magnitude of the new heat- 
ing term, describes the associated modifications to the thermal 
history of the IGM, and shows how this can explain some re- 
cent observations of the Lya forest. Pfrommer et al. (2012, 
hereafter Paper III) considers the impact the new heating term 
has upon the structure and statistics of galaxy clusters and 
groups, and upon the ages and properties of dwarf galaxies 
throughout the Universe, generally finding that blazar heat- 
ing can help explain outstanding questions in both cases. An 
additional follow-up paper by Puchwein et al. (2011) shows 
that when combined with most recent estimates of the evolv- 
ing photoionizing background and hydrodynamic simulations 
of cosmological structure formation, the heating by blazars re- 
sults in excellent quantitative agreement with observations of 
the mean transmission, one- and two-point statistics, and line 
width distribution of high-redshift Lya forest spectra. In par- 
ticular, these successes depends upon the peculiar properties of 
the blazar heating via the dissipation of plasma instabilities. 

For all of the calculations presented below we assume the 
WMAP7 cosmology with Iiq = 0.704, ^dm = 0.227, = 
0.0456, and Qa = 0.728 (Komatsu et al. 2011). 

2. THE FATE OF VERY HIGH ENERGY GAMMA RAYS AND THE 
PROPERTIES OF THE RESULTING PAIR BEAM 

2.1. Propagation and Absorption of Very High-Energy 
Gamma rays 

Sources of VHEGRs are necessarily attenuated by the pro- 
duction of pairs upon interacting with the EEL. Namely, when 
the energies of the gamma ray (£) and the EBL photon (/sebl) 
exceed the rest mass energy of the pair in the center of mass 
frame, i.e., 2£'£'ebl(1 -cos0) > 4-m^c'^, where 9 is the relative 
angle of propagation in the lab frame, an e^ pair can be pro- 
duced with Lorentz factor 7 ~ E /Im^c^ (Gould & Schreder 
1967). The optical depth the Universe presents to high-energy 
gamma rays depends solely upon the energy of the gamma ray 
and the evolving spectrum of the EBL'*. Detailed estimates for 
the EBL spectrum and its evolution have produced an estimated 
mean free path of TeV photons of 

^Pp(£,^) = 35(p^)'(i^)'Mpc, (1) 

where the redshift evolution is due to that of the EBL alone, is 
dependent predominately upon the star formation history, and 
C = 4.5 for z < 1 and C = for z > 1 (Kneiske et al. 2004; 
Neronov & Semikoz 2009)^. The resulting optical depth to 

As a consequence, careful studies of the absorption of high-energy gamma 
rays from extragalactic sources has been suggested as a way to probe the EBL 
(see, e.g., Gould & Schreder 1967; Stecker et al. 1992; Oh 2001). 

^ Despite the fact that the EBL contribution from starbursts peaks at z = 3 and 
declines rapidly afterward, galaxies and Type 1 active galactic nuclei compen- 




z 



Figure 1. Top; Pair-production optical depth of a gamma-ray photon emitted 
with energy E originating at redshift z as it propagates to the Earth, te . Bottom; 
Pair-production optical depth of a gamma ray of energy E propagating across a 
Hubble distance at redshift z,Tfj. In both cases the thick black line shows r of 
unity, blue and red lines show optically thick and thin regions of the parameter 
space, respectively. The sudden break in the contours at z = 1 is due to the 
assumed EBL evolution (see text). 

VHEGRs emitted with energy E from an object located at z is 

te(E,z)= f r ocE, (2) 

Jo D,,[E(l+z')/(l+z),z']Hiz')(l+z') 

(in which H{z) is the Hubble factor), is shown in the top panel 
of Figure 1.^ From this it is evident that above lOOGeV the 
Universe is optically thick to sources at z > 1 (cf. Frances- 
chini et al. 2008). A related, and perhaps more appropriate 
measure of the optical depth is that across a Hubble length, 
th = c/{DppH), shown in the bottom panel of Figure 1, pro- 
viding a sense of how opaque the Universe is as a function of 
redshift. In all cases it is clear that at E > lOOGeV photons 
will pair-produce on the EBL for z ^ 10. In the future this 
may not be the case, since for sufficiently small z, th decreases 
with decreasing z due to the dramatic decrease in the EBL pho- 
ton density associated with both the Hubble expansion and the 
slowing of star formation. 

2.2. Temperature of the Ultra-Relativistic Pair Beam 

The pairs produced by the TeV gamma rays are ultra- 
relativistic, with typical Lorentz factors of E /InieC^ ~ 
10*(£'/TeV). They also necessarily constitute a cold, highly 
anisotropic, dilute beam propagating through the IGM. This 
follows immediately from the intergalactic distances traversed 
by the gamma rays and the comparatively small EBL photon 
energies (as seen in the IGM frame). Here we estimate the 
properties of this plasma beam. 

sate for the lost flux until i = 1 . See, for example. Figure 3 from Franceschini 
et al. (2008). 

^ Note that the optical depth experienced by a photon that is absented to 
have energy E' is given by te[E'(\+z),z\. 
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Figure 2. Initial pair beam cooling rates due to the kinetic oblique instability 
(thick solid) and inverse Compton scattering (dotted) as a function of VHEGR 
energy (E) at a number of redshifts (z). In all cases 1+5 = 1, con'esponding to 
a mean-density region, and the isotropic-equivalent luminosity of the source at 
energy E, ELe, is 10'*^ ergs"', similar to the brightest TeV blazars seen from 
Earth. Finally, we list the initial pair Lorentz factor, 7, and cooling lengthscale 
along the top and right axes, respectively. 

The momentum dispersion of the resulting pairs is set by the 
gamma-ray spectrum, geometry of the TeV source, distribu- 
tion of the EBL photons, and heating due to pair production. 
Of these, only the last plays a significant role in setting the 
transverse momentum dispersion. The center-of-mass frame, 
i.e., the "beam frame", momentum dispersion resulting from 
pair production is roughly nieC. This results in an IGM-frame 

transverse momentum dispersion of p±/p\\ — 10"^ (£ /TeV) . 
With p±,meC <C p\\, the temperature'^ associated with this 

^ The contribution to the transverse momentum dispersion from the finite 
source size, due to the slightly different orientations of photons from opposing 
sides of the TeV emission region of size £, is p±/p\\ < £/2Dpp. With £ ~ 

10''*cm, implied by the X-ray variability timescales (Tramacere et al. 2009; 
Abdo et al. 2010d), this gives p±/p\\ < 5 X lO"'' [(1 +z)/2] (E/TeV)^. 
Similarly, since the creation of pairs is dominated by EBL photons near the 
pair-production threshold, the typical energy of the relevant EBL photons is 
roughly Am^c^/E (i.e., twice the threshold value for transverse EBL photons), 
and thus px/p|| ~ Amjc'^/E- ~ 1 X lO-'^ (£/TeV)"'. 

^ The temperature of an anisotropic particle distribution is inherently ill de- 
fined. Here we identify a temperature with the momentum dispersion of the 
beam using the drifting Maxwell- Juttner distribution employed by Bret et al. 
(2010a) (eq. 1 therein). 



transverse momentum dispersion is 

2 

-5x10" 

.P\\ ' 



'(tcv) 



(3) 



Since the transverse momentum dispersion of the pairs is much 
smaller than that associated with the bulk motion of the beam 
(i.e., since p±/p\\ ^ 1), we may safely assume that the beam 
is transversely kinematically cold. 

2.3. Density of the Ultra-Relativistic Pair Beam in General 

The density of the pair beam at a given point within the IGM 
is set by the rate at which pairs are produced, duration of the 
TeV emission, advection of the pairs through the IGM, and the 
processes by which they lose their kinetic energy. That is, the 
evolution of the density of pairs per unit Lorentz factor, is 
governed by the Boltzmann equation: 

c dr^n 



dr 



. dn^ 
97 



'7 ! 



(4) 



where the left-hand side assumes all the pairs are moving away 
from the TeV source relativistically (v'' = c and p'' = jniec), 
and the right-hand side corresponds to pair production. Gener- 
ally, we may neglect advection, which alters n^ over timescales 
of Dpp/c, much longer than any relevant timescale of interest 
here (i.e., {c / r^)dr~n^ / dr may be neglected). Furthermore, 
for most of the potential sources we will consider (primar- 
ily TeV blazars) we will assume that the duration of the TeV 
emission is sufficiently long that n^ reaches a steady state (i.e., 
dn^/dt = 0). In this case, we have '^dn^jd^ ~ ri^. 

Making further progress requires us to define the spectrum 
of the pairs, which itself depends upon the spectrum of the 
gamma rays and the energy dependence of the cooling pro- 
cesses. Nevertheless, we may obtain an estimate of the beam 
density in the vicinity of a given Lorentz factor, nt, ~ 7«^, by 
setting dn^/d^ ~ -n-y/7 = -ni,/^^. 

The source term is given by twice (since each gamma-ray 
produces two leptons) the rate at which high-energy gamma 
rays with energy E ~ lyngC^ annihilate within the region of 
interest, i.e., 7«'^ = 2{EdN/dE)/Dpp = IFg/Opp, where is 
the gamma-ray number flux, with units of photons cm~^s~'. 
Thus, upon defining a cooling rate F = -7/7, we have 



nt, ' 



2Fe 
Z)ppF' 



(5) 



i.e., the density of pairs of a given energy is determined by 
balancing cooling and pair creation. 

Generally, F is a function of energy and beam density, as 
well as external factors (e.g., seed photon density, IGM den- 
sity, etc.). Thus this gives a non-linear algebraic equation to 
solve for n/,, the particulars of which depend upon the various 
mechanisms responsible for extracting the bulk energy of the 
beam. In practice, given expressions for F, associated with the 
processes discussed in following section, we solve Equation (5) 
numerically to obtain ni,{E,FE,z). 

Which mechanism dominates the cooling of the beam de- 
pends upon a variety of environmental factors and the prop- 
erties of the pair beam itself. Nevertheless, inverse-Compton 
cooling via the cosmic microwave background (CMB) pro- 
vides a convenient lower limit upon F, and thus an upper limit 
upon ni,. This is a function of z and 7 alone, given by 



ic ■ 



4o-T«CMB 



-7: 



1.4x 10"2°(l+z)'*7 s"' 



(6) 
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where ctt denotes the Thompson cross section. The strong red- 
shift dependence arises from the rapid increase in the CMB 
energy density with z (mcmb fX (1 +zf ). Furthermore, since it 
is cx 7, inverse-Compton cooling is substantially more efficient 
at higher energies. The associated cooling rate is shown as a 
function of E for a number of redshifts in Figure 2. 

When we set F ~ Fic, we obtain the following upper limit 
upon the beam density: 



2Fe 



DppFic 27rZJ^pric 



3.7 X 10 



-22 



1+: 
2 



3C-4 



/ ELe 



V1045ergs-i J VTeV 



/ E 



(7) 

Where we have defined Le to be the isotropic-equivalent lu- 
minosity (per unit energy) of a source located a distance Dpp 
from the region in question. Setting ELe to a typical value 
(lO'^^ergs"') gives an idea of the typical pair-beam densities. 
Note that despite the large blazar luminosities we consider, the 
associated beams are exceedingly dilute, a point that is of crit- 
ical importance in the following section. Since Fic is indepen- 
dent of lib, this has no implication for inverse-Compton cooling 
itself. 

3. COOLING ULTRA-RELATIVISTIC PAIR BEAMS VIA PLASMA 
BEAM INSTABILITIES 

Plasma beams are notoriously unstable, with the instabili- 
ties driven by the anisotropy of the lepton distribution func- 
tion. Here we consider the implications of these instabilities 
upon the ultimate fate of the kinetic energy in the TeV-blazar 
driven pair beams. 

3.1. A Fundamental Limit Upon Plasma Cooling Rates 

For collective phenomena to be relevant, it is necessary for 
many pairs to be present within each wavelength of the unsta- 
bly growing modes. As we shall see, the relevant scale for the 
beam plasma instabilities we describe below is the plasma skin 
depth of the IGM, 



2t:c 

UJp 



I nm^c^ 
e^nicM 



= 2.3 X 10-^(1 +,5r'/\l +2)--'/- pc, (8) 



-3/2, 



where ujp = ^jA-Ke^nicu/me is the IGM plasma frequency and 
niGM = 2.2 X lQ~'^{l+5){l+zf cm~^ is the IGM free-electron 
number density. Generally, the growing mode must be uni- 
form on scales considerably larger than Ap, both longitudinally 
and transversely (otherwise it is not well-represented as a sin- 
gle Fourier component), and thus the volume in which many 
particles must be present is much larger than that defined by a 
sphere of diameter \p. Nevertheless, this gives us a conserva- 
tive constraint, i.e., we require 



-A^ni> 1 
6 



(9) 



With Equation (5) this gives a maximum plasma cooling rate: 



Fe 
-L plasma 



(10) 



Plasma processes with cooling rates that exceed this limit nec- 
essarily saturate near this cooling rate. Such a super-critical 
process can potentially operate only until the beam density is 




Figure 3. Limiting isotropic-equivalent luminosities as a function of redshift 
for various gamma-ray energies for 1 + (5 = 1 region of the IGM. Grey lines 
show the limit defining the applicability of the plasma prescription, below 
which less than a single beam lepton is found in a volume 7rAp/6. Below 
these collective phenomena are unimportant. Black lines show the luminosity 
at which linear plasma cooling rate begins to dominate inverse-Compton cool- 
ing. Both limits are sensitive functions of the IGM density, scaling as (1 +5)'/^ 
and (1+5)'/^, respectively. Thus in low/high density regions these limits be- 
come moderately more/less permissive. For reference, the sources listed in 
Table 1 are also plotted (at £ = 1 TeV), with HBL, IBL, radio galaxies, and 
quasar's shown by the the blue triangles, green squares, red hexagons and ma- 
genta circles, respectively. Filled points indicate sources that have been used 
to estimate the IGMF (see Section 4). The only objects that fail to meet the 
plasma-cooling luminosity criterion are the radio galaxies M87 and Cen A, 
both of which are detectable only due to their close proximity. 



driven below the value at which pairs can support collective 
phenomena, at which point they necessarily quench. How- 
ever, after plasma cooling ceases, the plasma beam density 
rises again (since Fic is always less than Fpiasma in practice), 
and thus the super-critical plasma cooling may resume. Since 
the efficiency of the plasma cooling decreases smoothly to zero 
at the critical density, this sequence stabilizes near Fpiasma- The 
associated excluded region is shown by the grey region in the 
upper-left corner of Figure 2. 

The upper limit upon nt obtained in Equation (7) implies an 
analogous limit upon the relevant source isotropic-equivalent 
luminosities. Again, this arises from requiring a sufficient 
number of beam pairs to be present to support plasma pro- 
cesses. Since the beam density is linearly dependent upon the 
source luminosity, this gives a constraint upon the latter: 

ELe»12E^'^^^'' 



= 3.3 X 10^'*(l-i-(5) 



l+Z 



8.5-3C 



TeV 



(11) 



effectively defining a luminosity cut-off, below which collec- 
tive plasma phenomena may be ignored. This limiting luminos- 
ity shown as a function redshift by the grey lines in Figure 3 for 
various gamma-ray energies. While the luminosity limit does 
depend upon z, it is clear that all of the observed TeV blazars 
(listed in Table 1) are sufficiently luminous for their resulting 
pair beams to support collective phenomena at the redshifts of 
interest for active galactic nuclei (AGNs). The only source to 
fall marginally below the limiting luminosity at 1 TeV is the ra- 
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dio galaxy Cen A (the lower-most point in Figure 3), the closest 
and dimmest object in Table 1. 

3.2. Cold Plasma Beam Instabilities 

Within astrophysical contexts there are at least two well 
known beam instabilities: two-stream and Weibel, the latter 
having been suggested as a mechanism for magnetizing strong 
shocks (Medvedev & Loeb 1999), and both implicated in the 
coupling at collisionless shocks (Spitkovsky 2008). These are, 
however, simply different limiting examples of the same under- 
lying filamentary instability, for which the maximum growth 
rate exceeds either - the so-called "oblique" mode, which we 
discuss below (Bret et al. 2004, 2005a; Lemoine & Pelletier 
2010). Note that since the pair beam is neutral, it contains its 
own return current and thus beam instabilities that arise due 
to the electron return currents within the background plasma 
(e.g., the Bell and Buneman instabilities, see Bret 2009) are 
not relevant. 

Here we discuss the nature of these instabilities, and their 
growth rates in the cold-plasma limit (i.e., mono-energetic 
beams). While the low-temperature approximation is unlikely 
to be applicable in practice for the beams of interest here, it 
provides a convenient limit in which to present the relevant pro- 
cesses within their broader context. For a similar reason, and 
because we have not found analogous derivations elsewhere in 
the astrophysical literature, we present the ultra-relativistic pair 
beam instability growth rates for the two-stream and Weibel in- 
stabilities in Appendixes A.l and A. 2, only summarizing the 
results here. We defer a discussion of the more directly rele- 
vant warm-plasma oblique instability to the following section. 
Generally, we find that the plasma instabilities are capable of 
dominating inverse Compton cooling as a means to dissipate 
the bulk kinetic energy of the pair beams from TeV blazars. 

The pair two-stream instability arises due to the interaction 
of the anisotropic electron and positron distribution functions 
with the comoving background electrostatic wave (i.e., k || p, 
where k is the electromagnetic mode wavevector and p is the 
beam momentum) with wavelength ~ Xp for n/,/«iGM ^ 1, 
which is generally the case of interest here. The associated 
cooling rate^ in the cold-plasma limit is 



me \2«igm/ 7' 



Tts ■ 



(12) 



which depends only weakly upon the IGM density and the pair 
beam density, though decreases rapidly as 7 becomes large. 

The Weibel instability is associated with coupling to a sec- 
ularly growing, anharmonic transverse magnetic perturbation 
(i.e., k _L p). The most rapidly growing wavelength is again 
that associated with the background plasma skin depth, ~ \p, 
with associated cooling rate in the cold-plasma limit 



Fw - 



\6iie^ni, 



(13) 



which depends only upon the pair beam density and the pair 
Lorentz factor At large 7 this is suppressed more weakly 
than the two-stream instability, though is a moderately stronger 
function of «;,. 

' Note that since we are interested in the rate at which energy is lost, these 
are necessarily twice the instability "growth" rates, defined by the e-folding 
time of the electromagnetic wave amplitudes. 



The computation of the growth rates of the two-stream and 
Weibel instabilities are greatly simplified by the particular ge- 
ometries of the coupled electromagnetic waves and beam mo- 
menta. However, a generalized oblique treatment has shown 
the presence of continuum of unstable modes (Bret et al. 2004, 
2005a; Bret 2009; Bret et al. 2010b; Lemoine & Pelletier 
2010), characterized by the orientation of their wave- vector 
relative to the bulk beam velocity. Of these neither the two- 
stream nor Weibel are generally the most unstable. Rather, 
generically the most robust and fastest growing mode occurs 
at oblique wave-vectors, and thus referred to as the oblique in- 
stability by Bret et al. (2010b). The cold-plasma cooling rate 
of this maximally-growing mode is 
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(14) 



This can be much larger than the two-stream and Weibel 
growth rates when «/, /higm ^ 1 and 7^1. 

3.3. Warm Plasma Beam Instabilities 

The cold-plasma oblique instability cooling rates dominate 
inverse-Compton cooling by orders of magnitude over the re- 
gion of interest. However, at the very dilute beam densities 
of relevance here, the cold-plasma approximation requires ex- 
ceedingly small beam temperatures. Above an IGM-frame 
temperature of roughly 
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pairs can traverse many wavelengths of the unstable modes 
over the cold-instability growth timescale, a situation com- 
monly referred to as the kinetic regime. As a consequence, 
significant phase mixing can occur, substantially reducing the 
effective growth rate (Bret et al. 2010a). 

The way in which finite beam temperatures limit the growth 
rate depends sensitively upon the nature of the velocity disper- 
sion and the modes of interest (see, e.g., Bret et al. 2005b). For 
example, the two-stream instability, associated with wave vec- 
tors parallel to the beam, enters the strongly suppressed "quasi- 
linear" regime when the parallel momentum dispersion is large, 
though is insensitive to even large transverse dispersions. Con- 
versely, the Weibel instability, associated with wave vectors 
orthogonal to the beam, is sensitive to even small transverse 
velocity dispersions but unaffected by large parallel velocity 
dispersions. This is simply because a given mode can tolerate 
large velocity dispersions within but not across the phase fronts 
of the unstable electromagnetic modes. For the situation of in- 
terest here, dilute beams and cool IGM (i.e., kTicM/nifC^ ^1), 
the oblique modes are nearly transverse, and thus sensitive pri- 
marily to large transverse velocity dispersions. 

Nevertheless, even for the small temperatures we have in- 
ferred for the pair beams, we find ourselves in the kinetic 
regime. In this case the oblique instability cooling rate has 
been numerically measured to be 
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where we have set the beam temperature in the beam frame to 
nieC^/k, and thus kTh = nifC^/'y (Bret et al. 2010a). Both the 
cold and hot growth rates have been verified explicitly using 
particle-in-cell (PIC) simulations, though for somewhat less di- 
lute beams than those we discuss here (Bret et al. 2010b). 

When the kinetic oblique mode dominates the beam cooling 
the beam density is given by 
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The associated cooling rate is then 
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This is a stronger function of gamma-ray energy than inverse- 
Compton cooling, implying that it will eventually dominate at 
sufficiently high energies, assuming a flat TeV spectrum. In 
addition it is a very weak function of 5, being only marginally 
faster in lower-density regions, and thus the cooling of the 
pairs is largely independent of the properties of the background 
IGM. 

The rates obtained by numerically solving Equation (5) for 
rtb, with r = Tic + rM,k, are shown for a number of redshifts in 
Figure 2. For the luminosity shown {ELe = 10'*'' ergs"', typical 
of the bright TeV blazars) at z < 4 plasma cooling dominates 
inverse Compton above a TeV. In the present epoch, Fm.Ic is 
roughly two orders of magnitude larger than Fjc for bright TeV 
blazars. 

The luminosity dependence of Fjvi.k implies a luminosity 
limit below which inverse Compton does dominate the linear 
evolution of the pair beam. 
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(note that at this luminosity F = Fic + Fm.Ic = 2FM,k, and thus rib 
is half the value shown in Equation (17)). This limit is shown 
as a function of redshift for a number of different energies 
by the black lines in Figure 3. Note that these lie above the 
corresponding limits associated with the applicability of the 
plasma prescription, suggesting that in practice the beam does 
support collective phenomena. The critical luminosity ranges 
from ^ lO'^'ergs"' to lO^^ergs"', depending upon redshift 
and energy of interest. It also depends upon the over-density 
as roughly cx (1 + 5)'/", and thus at low densities the critical 
luminosity is moderately smaller. The only two sources in 
Table 1 which fall below the plasma-cooling luminosity limit 
are the radio galaxies M87 and Cen A, both of which are 
detectable only as a result of their close proximity. At 1 TeV, 
all but a handful of the remaining sources lie more than an 
order of magnitude above this limit, and in the case of the two 
sources that dominate the TeV flux at Earth, more than two 
orders of magnitude above this limit. 



3.4. Intuitive Picture of the two-stream, Weibel, and Oblique 
Instabilities 

Given the technical nature of our discussion above, it is use- 
ful to have a qualitative understanding of these instabilities. We 
caution, however, that intuitive pictures of plasma processes 
frequently fail to capture all the relevant physics. Hence, gen- 
eralization of these intuitive pictures beyond their limited range 
of applicability is potentially misleading. This is explicitly il- 
lustrated by our examples here: all of the instabilities discussed 
in this work belong to the same family, i.e., instabilities that 
arise from interpenetrating plasmas, but the underlying quali- 
tative pictures for each differ substantially. 

We begin with the Weibel (or filamentation) instability 
(Weibel 1959), for which a mechanical viewpoint (an initial 
magnetic perturbation deflects particles into opposing current 
streams that reinforce the perturbed field) can be found in 
Medvedev & Loeb (1999), to which we refer the interested 
reader. Here we present picture that although having the virtue 
of being simpler is not entirely correct: because like cur- 
rents attract, small-scale current perturbations arising out of 
the fluctuations within the interpenetrating plasmas will coa- 
lesce preferentially to produce increasingly larger-scale cur- 
rents. These induce stronger magnetic fields, and thus larger 
attractive Lorentz forces between neighboring currents; a pos- 
itive feedback loop develops leading to instability. This pro- 
cess continues until the associated magnetic field strengths be- 
come sufficiently large to disrupt the currents (in the case of 
equal density beams) or until the transverse velocity of the con- 
stituent particles is large enough to efficiently migrate between 
current structures on the linear growth timescale, i.e., enter the 
kinetic regime. 

We should note that the aggregation of currents in the Weibel 
instability does not rely upon oscillatory waves. Hence, there 
is no oscillatory component to this instability - it is a purely 
growing mode. In contrast, the two-stream instability is an 
overstable mode, where the oscillatory components are Lang- 
muir (or plasma) waves with a wavevector parallel to the beam 
velocity. These are longitudinal waves, associated with local 
charge oscillations, and are completely described by a propa- 
gating perturbation in the electrical potential. As particles in 
the beam traverse Langmuir waves in the background plasma, 
they experience successive periods of acceleration and decel- 
eration, with electrons (positrons) collecting in minima (max- 
ima) of the electric potential, where the particle speeds are at 
their smallest. ' ' This charge-separated bunching of the beam 
plasma enhances the background electric perturbation, poten- 
tially growing the background Langmuir wave. 

The bunching within the beam is simply an excitation of 
Langmuir waves within the beam plasma itself. Thus the 
growth of the charge perturbations in the background and beam 
plasmas corresponds to the resonant coupling between Lang- 
muir waves in the background and beam. When the beam den- 
sity is much less than the background density, as is the case 
here, background and beam Langmuir waves only overlap in 
frequency, and therefore satisfy the conditions for resonance, 
when the wavevector of the latter is parallel to the beam veloc- 

The linear analysis of the Weibel instability assumes wave-like distur- 
bances in the plasma. However, there is no associated restoring force to these 
waves and, hence, they do not oscillate. 

' ' One tempting aspect is to think of these particles as almost massless and 
hence imagining that these electric fields are quickly shorted out. While this is 
the case in most astrophysical process, we urge the reader to resist this tempta- 
tion because these instabilities occur on short enough timescales that the mass 
of these charged particles plays a crucial role in the physics of the instability. 
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ity (see the discussion above Equation (A9)). This is always 
satisfied for the family of comoving beam Langmuir waves 
(i.e., waves which in the beam frame move counter to the back- 
ground plasma). However, of particular importance for the 
two-stream instability are the counter-propagating beam Lang- 
muir waves (i.e., waves which in the beam frame move parallel 
to the background plasma). If the beam velocity exceeds the 
phase velocity of these waves (as seen in the beam frame), the 
counter-propagating wave will be dragged in the direction of 
the beam (as seen in the background frame), and therefore has 
a wavevector which satisfies the resonant condition. 

Nevertheless, the counter-propagating Langmuir wave still 
carries momentum in the direction opposite to the beam. Thus, 
as the counter-propagating wave grows, the momentum, and 
therefore energy, of the beam-wave system necessarily de- 
creases. This implies that the counter-propagating Langmuir 
wave is also a negative energy mode as seen in the background 
frame. As a consequence, the resonant coupling can trans- 
fer energy to the positive-energy background wave from the 
negative-energy beam wave, while growing the amplitudes of 
both, and thereby leading to instability. Note that if the phase 
velocity of the counter-propagating Langmuir wave is larger 
than the beam velocity, it no longer is dragged in the direc- 
tion of the beam and no longer satisfies the necessary resonant 
condition. For a distribution of particles, this constraint upon 
the velocities within the beam corresponds to the familiar Pen- 
rose criterion (Sturrock 1994; Boyd & Sanderson 2003), and is 
satisfied in the pair beams resulting from VHEGRs. 

The oblique instability encompasses the two-stream and 
Weibel instabilities, though the most unstable mode is most 
similar to the former in that the intuitive picture focus solely 
on the electrostatic forces, ignoring electromagnetic forces.'^ 
In practice, this is a relatively good approximation and can be 
used to calculate the growth of these modes in idealized sit- 
uations (e.g. Nakar et al. 2011; Bret et al. 2004). The qual- 
itative picture proceeds similarly to that for the two-stream 
instability described above, with the minor modification that 
the perturbing background Langmuir waves now move at an 
angle 6 relative to the relativistic beam. As a consequence, 
the resonant Langmuir waves have a phase velocity such that 
Vk ~ ccos6'. From the simple intuitive picture above, if 9 is 
selected such that the projected beam velocity is slightly faster 
than the nearly resonant Langmuir wave, an instability devel- 
ops. 

Finally, to understand why the growth rates between the two- 
stream instability above and the oblique instability differ, it is 
useful to make the approximation that the electric fields gener- 
ated are in the direction of the k-vector In the two-stream case, 
the electric field must slow down (or speed up) particles. This 
gets progressively harder for more relativistic particles. In the 
oblique case, the electric field deflects the particles, changing 
their projected velocity. While this is also more difficult for 
more relativistic particles, this is not nearly as hard as chang- 
ing the particles parallel (along the beam) velocity. Hence, the 
oblique instability more easily drives charge density enhance- 
ments (and therefore instabilities) at large 6, i.e., easier deflec- 
tion, than the two-stream instability. 

3.5. Non-Linear Saturation 

We have thus far only treated the linear development of the 
relativistic two-stream, Weibel, and oblique instabilities. How- 

The picture we present here for the obhque instabihty is discussed in 
Nakar etal. (2011). 



ever, the impact pair beams have upon the IGM, gamma-ray 
cascade emission, and measures of the IGMF will ultimately 
depend upon their nonlinear development. To address this, 
however, we are presently forced to appeal to analytical and 
numerical calculations of systems in somewhat different (and 
less extreme) parameter regimes. 

Motivated by the applicability of the Weibel instability in 
the context of GRBs, the nonlinear saturation of the relativis- 
tic Weibel instability for equal density plasma beams is well 
understood analytically and numerically. Initially, the Weibel 
instability rapidly grows until the energy density of the gener- 
ated magnetic field becomes of order the kinetic energy of the 
two beams (Silva et al. 2003; Frederiksen et al. 2004; Chang 
et al. 2008). Analytically, Davidson et al. (1972) argued that 
the Weibel instability would saturate when the generated mag- 
netic fields become so large that the Larmor radius of the beam 
particles is of order the skin depth, i.e., when the energy of 
generated magnetic fields is equal to the kinetic energy of two 
equal density relativistic beams (see also Medvedev & Loeb 
1999). The particles rapidly isotropize with a Maxwellian dis- 
tribution (Spitkovsky 2008), i.e., heat, and the magnetic en- 
ergy then rapidly decays within an order of a few tens of skin 
depths (Chang et al. 2008). Hence for two equal density, rela- 
tivistic, interpenetrating beams, the Weibel instability converts 
anisotropic kinetic energy into heat. However, as we have al- 
ready noted, for the pair beams of interest here the Weibel 
instability is completely suppressed for tiny transverse beam 
temperatures. Hence, while this instability may initially oper- 
ate, it may quickly become suppressed if it results in significant 
transverse heating of the beam. 

Unlike the Weibel instability, the two-stream and kinetic 
oblique instabilities continue to operate, though more slowly 
than implied by their cold-plasma limits, in the presence of sub- 
stantial beam temperatures. Due to the geometry of the coupled 
modes, the oblique instability is primarily sensitive to trans- 
verse beam-velocity dispersions, though shares the resistance 
to beam heating with its two-stream cousin (Bret et al. 2005b; 
Bret 2009; Bret et al. 2010b; Lemoine & Pelletier 2010). Un- 
like the two-stream instability, the oblique instability in the pa- 
rameter range of interest here is also largely insensitive to lon- 
gitudinal heating. Thus, generally, it appears that the oblique 
instability is substantially more robust than its more commonly 
discussed brethren. 

What is less clear a priori is if these instabilities primarily 
heat the beam or primarily heat the background plasma. Here 
we appeal to the numerical simulations of Bret et al. (2010b), 
where a mildly relativistic beam (7 = 3) penetrating into a hot, 
dense background plasma (beam-to-background density ratio 
of 0.1) was studied. In these the oblique instability resulted 
in a significant fraction (~ 20%) of the beam energy heating 
the background plasma before the heating of the beam sup- 
pressed the oblique instability in favor of the two-stream insta- 
bility. The relative effectiveness with which the beam heats the 
background plasma is due to the efficiency with which the lon- 
gitudinal electrostatic modes are dissipated in the background 
plasma; unlike electromagnetic modes (e.g., those generated 
by the Weibel instability), electrostatic modes are rapidly dis- 
sipated via Landau damping. 

We note that Bret et al. (2010b) found that the heating of the 
beam by the oblique instability eventually led to its suppres- 
sion, allowing the two-stream instability to grow, continuing 
the dissipation of the beam kinetic energy. As a result, in their 
simulations a total of 30% of the energy was deposited into the 
background plasma via a combination of the oblique and two- 
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stream instabilities, i.e., an additional 10% of the beam energy 
was thermalized via the two-stream instability during and after 
the suppression of the oblique instability. 

In our case, we expect that much more beam energy (more 
than the 20%, and possibly up to ~ 100%) will be deposited 
into the background IGM because we are much deeper into the 
regime in which the oblique instabiUty dominates, i.e., 7^1 
and rib/niQM <C 1. In particular, for the case of interest here, 
7 - lO*" and 
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The extremely large Lorentz factor and tiny density ratio make 
it computationally prohibitive to assess the beam evolution 
with numerical PIC directly. Nevertheless, because we find 
ourselves in a regime in which the oblique instability is much 
more strongly dominant than that simulated in Bret et al. 
(2010b), we expect the linear growth of the kinetic oblique in- 
stability to continue for much longer before the beam changes 
character and moves out of the oblique-dominated regime — 
potentially once ~ 100% of the beam energy has been dissi- 
pated. However, this remains to be studied in future work. 

The effect of nonlinear processes might also effect the evo- 
lution of the linear instability. For instance, Lesch & Schlick- 
eiser (1987) argued that for the relativistic electrostatic two- 
stream instability, nonlinear coupling to daughter modes arrest 
the growth of the linearly unstable mode at a very low mode 
energy. Hence, they claim that the electrostatic two-stream in- 
stability can only bleed energy from the beam at a slow rate. 

Utilizing the order of magnitude estimates in Lesch & 
Schlickeiser (1987) or the expression for nonlinear Landau 
damping in Melrose (1980), we have found that damping of 
the pair beam via the relativistic two-stream instability could 
be highly suppressed. The oblique instability may be similarly 
suppressed, however, this effect is much more marginal due 
to the instabilities much larger growth rate. Nevertheless, de- 
termining its behavior for the parameters relevant here is an 
important unanswered question that is left for future work. For 
the purposes of this paper, we will rely upon the intuition de- 
veloped from numerical studies of the oblique instability and 
argue that the beam ends up heating the IGM primarily. 

In what follows, we will presume the nonlinear evolution of 
the "oblique" or related plasma instabilities lead to the heating 
of the background IGM, i.e., beam cooling. However, we note 
that the beam itself may be the primary recipient of this kinetic 
energy, i.e., beam disruption. Most of our results depend 
critically on beam cooling and not beam disruption. This 
includes the effect of blazar heating on the IGM temperature- 
density relation studied in Paper II, the effects on structure 
formation studied in Paper III, the excellent reproduction of 
the statistical properties of the high-redshift Lya forest found 
in Puchwein et al. (2011), and the implications on the EGRB 
and the redshift evolution of TeV blazars studied in Section 
5. However, our conclusions on the inapplicability of IGMF 
constraints determined from the non-observation of GeV 
emission from blazars still remains in the presence of beam 
disruption. This is because the self-scattering of pairs in the 
beam would suppress this GeV emission in similar manner to 
an IGMF 



3.6. Suppression of the Oblique Instability by an IGMF 

The beam instabilities we have discussed have been analyzed 
primarily within the context of unmagnetized plasmas. How- 
ever, for a variety of theoretical reasons a weak IGMF is not 
unexpected. For example, a field strength of 10"'^ G is suf- 
ficient to explain the observed nG galactic fields via compres- 
sion and winding alone. Here we consider the implications that 
an IGMF has for the instability growth rates we have described 
above. 

A strong IGMF causes the ultra-relativistic pairs to gyrate, 
and therefore to isotropize, suppressing the growth of insta- 
bilities that feed upon the beam anisotropy (e.g., those we 
have described above). However, for this to efficiently quench 
the growth of the plasma beam instabilities, this isotropiza- 
tion must occur on a timescale comparable to the instability 
growth time, i.e., the Larmor frequency must be comparable 
to the cooling rate, F. This condition gives a lower-limit upon 
IGMF strengths sufficient to appreciably suppress the growth 
of plasma beam instabilities of 
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considerably larger than those typical of both, primordial for- 
mation mechanisms (Widrow 2002) and implied by galactic 
magnetic field estimates, assuming galactic fields are produced 
by contraction and winding alone. 

Because the VHEGRs emitted by the TeV blazars travel cos- 
mological distances prior to producing pairs, an IGMF capable 
of suppressing the plasma beam instabilities must necessarily 
have a volume filling fraction close to unity. In particular, it 
must permeate the low-density regions, where most of the cool- 
ing occurs. However, the Alfven velocity within those areas is 
extraordinarily small, roughly 
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is considerably larger, implying convection is much more effi- 
cient. Nevertheless, even after substantial heating via the ther- 
malization of the TeV blazar emission (Paper II) is taken into 
account, magnetic fields will have propagated < lOOkpc over a 
Hubble time via convection and much less via diffusion, imply- 
ing that a pervasive, sufficiently strong magnetic field can not 
be produced via ejection from galactic dynamos. While galac- 
tic winds can produce much faster outflows, vw ~ lO^'km s"', 
unless they inject a mass comparable to that contained in 
the low-density regions over a Hubble time, they are rapidly 
slowed via dissipation at shocks in the IGM, again limiting 
the spread of galactic fields. Moreover, we note that since the 
magnetic field must be volume filling to suppress the plasma 
beam cooling substantially it is insufficient to produce pockets 
of strong fields, and thus any galactic origin powered by winds 
requires a nearly complete reprocessing of the low-density re- 
gions. However, the ultimate thermalization of such fast, dense 
winds would raise the IGM temperature to ~ 10** K, in conflict 
with the Lya forest data and exceeding the entire bolometric 
output of quasars by at least a factor of two. For these reasons 
we conclude that a volume-filling strong IGMF would demand 
a primordial origin. 
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4. IMPLICATIONS FOR IGM MAGNETIC FIELD ESTIMATES 

The existence of plasma processes that can cool the pair 
beams associated with TeV blazars has profound consequences 
for efforts to constrain the IGMF using the spectra of TeV 
blazars. Here we describe how the reported IGMF limits have 
been obtained, the consequences of plasma cooling for these, 
and potential strategies for overcoming the constraints it im- 
poses. 

The general argument made in efforts to constrain the IGMF 
based upon the GeV emission from blazars proceeds as fol- 
lows: The beamed TeV blazar emission pair-creates off of the 
EBL. The resulting pairs subsequently up-scatter CMB pho- 
tons to GeV energies. In principle, this should produce an 
observable GeV excess, or bump, in the spectra of these ob- 
jects. However, in the presence of a large-scale IGMF, the 
ultra-relativistic pairs can be deflected significantly, directing 
the beamed GeV emission away from Earth. Thus, it is argued, 
the lack of a discernible GeV bump in a number of TeV blazars 
implies a lower limit upon the IGM field strength (as a func- 
tion of TeV jet opening angle and variability timescale). The 
crucial components of the argument are 

1. The TeV emission is beamed with the typical opening 
angles inferred from radio observations of AGN jets. 

2. The variability timescale within the TeV source is long 
in comparison to the geometric time delays between 
the original TeV and inverse-Compton produced GeV 
gamma rays due to the orbit of the pairs through some 
angle ?9, Af - 10^(Dpp/80Mpc)('(?/0.1 rad)2yr (see Der- 
mer et al. 2011). 

3. The pairs produced by TeV absorption on the EBL cool 
primarily via inverse-Compton scattering the CMB (and 
therefore evolve only due to inverse-Compton cooling 
and orbiting within the large-scale IGMF). 

The first of these is supported indirectly by the lack of TeV 
emission from non-blazars. The implied TeV source stability 
required by the second is at odds with the variability observed 
in blazars at longer wavelengths. However, since the TeV-GeV 
delay is a strong function of deflection angle, given an empir- 
ical limit upon the TeV blazar variability timescale (presently 
4yr, Dermer et al. 2011) it is possible to produce a substan- 
tially weaker constraint upon the IGMF (~ 10"'^ G) (Dermer 
et al. 2011; Taylor et al. 201 1; Takahashi et al. 2012). 

Plasma cooling provides a fundamental limitation for these 
methods, however, by violating the third condition explic- 
itly. Figures 3 and 4 imply that for all but the dimmest and 
highest-redshift (z > 4) gamma-ray blazars only a small frac- 
tion of the pair energy is lost to inverse-Compton on the CMB 

(/ic = Tic/ (Fic + rivi.k))-''' In particular, the black lines in 
Figure 3 shows where /ic = 0.5, i.e., roughly 50% of the TeV- 
photon power is ultimately converted into heat via plasma in- 
stabilities. As a result, the putative GeV component is typically 
much less luminous than otherwise expected, reducing the sig- 
nificance of non-detections substantially. 

This may be seen explicitly in Figure 4, which shows fic as 
a function of gamma-ray flux for E = 1.85 TeV (correspond- 
ing to a Comptonized-CMB photon energy of approximately 
3 GeV), at number of source redshifts. Typical values for the 

Here we implicitly assume that the nonlinear saturation of the plasma 
instabilities extract energy from the pair beam at the linear growth rate. 
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Figure 4. Fraction of the energy in the pair beam lost via inverse- 
Comptonization of the CMB at a number of different redshifts. In all cases 
the injection photon energy (as would have been observed at z = 0) is £ = 
1.85 TeV, resulting in 3 GeV up-scattered CMB photons, the energy at which 
the Fermi/LAI instrument is most sensitive. For reference, the sources listed 
in Table 1 are also plotted (at E = \ TeV), with HBL, IBL, radio galaxies, and 
quasars shown by the the blue triangles, green squares, red hexagons and ma- 
genta circles, respectively. Filled points indicate sources that have been used to 
estimate the IGMF, coiTesponding to (in increasing flux) PKS 0548-322, lES 
0347-121, lES 1101-232, RGB J0152+017, lES 0229+200, lES 1218+304, 
RGB J07 10+591, and Mkn 501. For reference, the inferred isotropic luminos- 
ity is shown in the top axis for sources located at z = 0.1. 

TeV blazars collected in Table 1 (including those employed by 
Neronov & Vovk 2010; Tavecchio et al. 2010, 2011; Dermer 
et al. 2011; Taylor et al. 2011; Takahashi et al. 2012; Dolag 
et al. 2011) lie in the range lO""* to 3 x 10"-, implying corre- 
spondingly small GeV Comptonization signals. 

More importantly, when the beam evolution is dominated 
by plasma instabilities, over the inverse-Compton cooling 
timescale the pair distribution necessarily becomes isotropized. 
As a consequence, the angular distribution of the resulting GeV 
gamma rays (i.e., the orientation of the GeV "beam") is no 
longer indicative of the beam propagation through a large-scale 
magnetic field. That is, one expects large-angle deviations re- 
gardless of the IGMF strength. 

Unfortunately, subject to the caveats of the preceding sec- 
tion, plasma instabilities appear to be the dominant cooling 
mechanism for the pair beams associated with the TeV blazars 
that have been used to constrain the IGMF thus far The 
isotropic-equivalent luminosities of these sources range from 
2.5 X 10"*"' ergs"' to lO'^^ergs"', placing them well within the 
plasma-instability dominated regime. As a result, the reported 
IGMF limits inferred from TeV blazars are presently unreli- 
able. 

Nevertheless, it may be possible to avoid the limitations im- 
posed by plasma instabilities. At sufficiently low pair densi- 
ties the cooling rates associated with the beam instabilities de- 
scribed in Section 3 fall below that due to inverse Compton. 
Moreover, at some point the plasma prescription breaks down 
altogether, suggesting that any plasma instabilities that operate 
on the skin-depth of the IGM are strongly suppressed. There 
are two distinct ways in which low pair densities can arise; low 
intrinsic luminosities and very short timescale events. 

Below isotropic-equivalent luminosities of roughly ~ 
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10 ergs" inverse-Compton cooling dominates the beam in- 
stabilities we describe in Section 3 at ~ 1 TeV. Below 
lO^'ergs"^ the plasma prescription itself breaks down alto- 
gether. While no TeV blazars with isotropic-luminosities be- 
low 10^*^ ergs"' are known, there are a handful of very nearby 
sources below lO'*" erg s"' . The two dimmest sources, the radio 
galaxies M87 and Cen A, however, are observable only due to 
their proximity (i.e., they have proper distances ^ £>pp), intrin- 
sically preventing them from providing a significant constraint 
on the IMGF. The highest source in Figure 4, and thus presum- 
ably the source with the largest fractional inverse-Compton sig- 
nal, is PKS 2005-489, a dim, soft TeV blazar at z = 0.071 . Even 
in this case, however, the inverse-Compton cooling timescale 
is roughly 4 times longer than the plasma cooling timescale. 
Hence, probing the IGMF with TeV blazars will require ob- 
serving considerably dimmer objects than those used thus far. 
Since doing so will likely require substantial increases in de- 
tector sensitivities'"^ we will not consider this possibility any 
further here. 

Low pair densities may also be produced by Umiting the du- 
ration of the VHEGR emission. Equation (5) was derived as- 
suming that the pairs had reached a steady state between their 
formation via VHEGRs annihilating upon the EBL and their 
cooling via inverse-Compton and plasma processes. However, 
it takes roughly a cooling timescale for the pair beam densities 
to saturate at this level, which can be as long as 10^ yr at some 
E and z- Therefore, generally there is a lag between the onset of 
the VHEGR emission and the time at which plasma processes 
begin to dominate the cooling of the beam. 

We can estimate this by setting nt, ~ 2f£Af/£)pp where At is 
the duration of the emission. This is the case when the pair den- 
sity is initially zero (i.e., it has been many cooling timescales 
since the last period of VHEGR emission) and F Af ^ 1 . In- 
serting this into the various cooling rates and setting rM,k = Fic 
gives an estimate for the maximum duration for which the pair 
density remains sufficiently low that inverse-Compton cooling 
dominates the cooling: 

Ar<5.3(l + 5y/^(^l±^j 

which may be found for the observed TeV sources in Table 1 . 
Typically Af ~ 300 yr for the blazars that have been used to 
constrain the IGM, though Af ranges from 10^ yr to 10"^ yr for 
TeV blazars generally. Note that it is not sufficient for emission 
to vary on Af; such variations will be temporally smoothed 
on the much larger cooling timescale. Rather, it is necessary 
for the source to have been quiescent (i.e. isotropic-equivalent 
luminosity considerably less than 10^^ ergs"') for periods long 
in comparison to F"' and have turned on less than a time Af <C 
F"' ago. As a consequence, only IGMF estimates using new 
TeV blazars can potentially avoid plasma instabilities in this 
way. 

A natural example of a class of transient gamma-ray 
sources is gamma-ray bursts (GRBs). For GRBs an isotropic- 

''^ As described in Section 5.1.2, the incompleteness of tlie TeV observa- 
tions are well-described by a single flux limit, corresponding to an isotropic- 
equivalent luminosity limit at 10-Mpc of roughly 10'*' ergs"'. The flux Hmit 
of Fenni is estimated to be roughly 4 X 10'*^^ ergs"' at z ~ 0. 1 (see Figure 23 
of Abdo et al. 2010c). 



equivalent energy of roughly 10 erg is required for the plasma 
instabilities to grow efficiently, comparable to the total ener- 
getic output of the brightest events. Using GRBs to probe the 
IGMF has been suggested previously, and limits based upon 
the lack of a delayed GeV component in some GRBs already 
exist, finding field strengths above ^ 10"^^ G (see, e.g., Plaga 
1995; Dai & Lu 2002; Guetta & Granot 2003; Takahashi et al. 
2008). Future efforts should be able to probe field strengths 
of 10"'5G for z < 0.2 (see Figure 4 of Takahashi et al. 2008). 
This is predicated, however, upon the existence of a signifi- 
cant VHEGR component in the prompt or afterglow emission 
of GRBs. 

Ando & Kusenko (2010) have attempted to directly detect 
the GeV signal in the vicinity of AGNs that do not exhibit 
TeV emission. That is, since an IGMF can in principle induce 
large-angle deviations in the propagation direction of the up- 
scattered GeV photon relative to the original VHEGR, one may 
search for a GeV excess at large distances from non-blazars. 
This is made much more difficult by the extraordinarily low 
surface gamma-ray brightness due to the large Dpp, random 
source orientations, and potential dilution associated with the 
beam diffusion. Nevertheless, Ando & Kusenko (2010) claim a 
marginal detection of this component in stacked Fermi images, 
though the Fenni point-spread function appears to be capable 
of completely explaining their result (Neronov et al. 2011). 
If such a GeV-halo were found around TeV-bright sources 
with the anticipated luminosity, it would argue strongly against 
plasma instabilities dominating the beam evolution. However, 
we note that for many Te V-dim sources inverse-Compton cool- 
ing still dominates the beam evolution, and thus in stacked im- 
ages of these sources GeV-halos may still be present. 

Note that, as discussed in Section 3.6, a strong IGMF (^ 
10"'^ G) could, in principle, suppress the growth of the plasma 
beam instabilities significantly. As we argue there, a suffi- 
ciently pervasive and strong IGMF would necessarily be pri- 
mordial in origin. However, few constraints upon such a pri- 
mordial field exist, especially since we have argued that the 
lack of an inverse Compton GeV bump cannot be used as a 
constraint on the IGMF. Hence, we are left with constraints 
from Big Bang Nucleosynthesis (B < 10"^ G, Kandus et al. 
2011, and references therein), from the CMB (B < 10""^ G, 
Barrow et al. 1997), and limits on the Faraday rotation of 
distance radio sources (B < 10"'' G, Vallee 2011, and refer- 
ences therein). A stronger limit can be obtained from limits 
upon the deflection angle (< 2°) of ultra-high energy cosmic 
rays {E > lO^^eV) from distant sources at a distance d > \b, 
which upon assuming a magnetic correlation length, A^, gives 

B < 10"" (£/1020eV) (Afi/lOOMpc)"' (tZ/A/;)"'^^ G (see for 
instance Waxman & Miralda-Escude 1996; de Angelis et al. 
2008), though this depends strongly upon the assumed struc- 
ture of the IGMF and small-scale fields may be significantly 
stronger than this (Kotera & Lemoine 2008). Proposed mech- 
anisms by which a primordial IGMF can be theoretically gen- 
erated are essentially unconstrained, yielding B-field strengths 
between B ~ 10"'-10"^''G (Kandus et al. 201 1). 

However, in the remainder of this work and in Paper II and 
Paper III, we show that the effect of the "oblique" (or simi- 
lar) instability may potentially explain many different observed 
phenomena in cosmology and high energy astrophysics. Tak- 
ing the beam instabilities we present here at face value, that so 
many observational puzzles can be simultaneously explained 
by the local dissipation of TeV blazar emission would im- 
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ply the strongest upper limit upon primordial magnetism to 
date. Namely, the apparent impact of TeV blazars upon the 
large-scale cosmological environment places a constraint on 
the IGMF ofB < lO-'^G. 

5. IMPLICATIONS FOR THE BLAZAR LUMINOSITY FUNCTION 

The EGRB corresponds to the diffuse gamma-ray back- 
ground (E > lOMeV) not associated with Galactic sources. 
This has been measured at high energies by EGRET (Sreeku- 
mar et al. 1998; Strong et al. 2004) and more recently strongly 
constrained between 200 Me V and lOOGeV by Fermi (Abdo 
et al. 2010b, a). The most recent Fermi estimate of the EGRB 
is a featureless power-law, somewhat softer than that found by 
EGRET, and does not exhibit the localized high-energy excess 
claimed by a reanalysis of the EGRET data (Strong et al. 2004). 

While bright, nearby blazars are resolved, and therefore 
excluded from the EGRB, the remainder is thought to arise 
nearly exclusively from distant gamma-ray emitting blazars, 
upon which the Fermi EGRB places severe constraints (Abdo 
et al. 2010b). Beyond z — 0.25 Fermi cannot detect blazars 
with isotropic-equivalent luminosities comparable to the ob- 
jects listed in Table 1, and thus the vast majority of such objects 
contribute to the Fermi EGRB. 

A significant EGRB above lOOGeV is not expected due 
to annihilation upon the EBL. However, if they operate effi- 
ciently, ICCs can reprocesses the VHEGR emission of distant 
sources into the Fermi EGRB energy range (i.e., < lOOGeV). 
Thus, a number of efforts to constrain the VHEGR emis- 
sion from extragalactic sources based upon the EGRB can be 
found in the literature (Narumoto & Totani 2006; Kneiske & 
Mannheim 2008; Inoue & Totani 2009). These have typically 
found that the comoving number density of VHEGR-emitting 
blazars could not have been much higher at high-z than it is to- 
day. Even with moderate evolutions (e.g., that in Narumoto & 
Totani 2006; Inoue & Totani 2009) consistent with the EGRET 
EGRB, substantially over-produce the Fermi EGRB, and are 
therefore believed to be excluded (Venters 2010). 

However, here we show that this conclusion is predicated 
upon the high-efficiency of the ICC. We have already shown 
that for bright VHEGR sources plasma beam instabilities ex- 
tract the kinetic energy of the first generation of pairs much 
more rapidly than inverse-Compton scattering. As a conse- 
quence, the ICC is typically quenched, substantially limiting 
the contributions of these sources to the EGRB. Thus, even 
with a dramatically evolving blazar population, e.g., similar to 
that of quasars, it is possible for these objects to be consis- 
tent with the Fermi EGRB. While our discussion of the EGRB 
and the evolution of blazars is predicated on the extraction of 
the pair-beam kinetic energy by plasma beam instabilities, our 
conclusions are not. Rather they will continue to hold as long 
as the pair beams are locally dissipated - via plasma beam in- 
stabilities or some other equally powerful mechanism. 

5.1. An Empirical Estimate of the TeV Blazar Luminosity 
Function 

Blazars dominate the extragalactic gamma-ray sky, and thus 
represent the best studied VHEGR source class. Since we are 
interested exclusively in those objects responsible for the bulk 
of the VHEGR emission, we necessarily concentrate upon the 
subset of blazars that are luminous VHEGR emitters. Of the 
28 TeV sources listed in Table 1, 23 are peaked at very-high 
energies (the HBL and hard IBL sources, for a full definition 
see below), and comprise what we will call collectively TeV 



blazars. This necessarily is limited to z ~ 0.1 as a consequence 
of the annihilation of VHEGRs upon the EBL. Nevertheless, 
we will find that this is remarkably similar to the local quasar 
luminosity function {(pg), and thus attempt to construct a blazar 
luminosity function by analogy: 



iizt/log[()L 



(25) 



in which is the comoving number density of blazars at red- 
shifts less than z with isotropic-equivalent luminosities below 
L. This is different from previous efforts to empirically con- 
strain 4>B from the EGRB (see, e.g., Narumoto & Totani 2006; 
Inoue & Totani 2009) in at least two ways. First, we begin with 
an empirically determined local 0b and attempt to extend this 
by placing the TeV blazars within the broader context of accret- 
ing supermassive black holes, rather than beginning with the 
EGRB and working backwards to infer an acceptable (j>B- Sec- 
ond, we are primarily concerned with (pB of TeV blazars, and 
specifically do not consider the contributions from other kinds 
of objects. While this does not represent a significant oversight 
in terms of the high-energy contributions to the EGRB, which 
almost certainly arises from the VHEGR-emitting blazars, it 
does mean that our conclusions regarding the luminosity func- 
tion of the TeV blazars do not necessarily apply to all Fermi 
AGN (e.g., the FSRQs, see below). 

5.1.1. Placing TeV blazars in Context 

The TeV blazars presumably fit within the broader context 
of the blazars observed by Fermi specifically, and AGNs gen- 
erally. For this reason, here we briefly review the physical clas- 
sification scheme based on the widely accepted AGN standard 
paradigm that provides a unified picture of the emission emis- 
sion properties of these objects (e.g., Urry & Padovani 1995). 
Specifically, we summarize the classes of objects believed to 
be capable of producing significant TeV luminosities and the 
potential physical processes responsible for the observed emis- 
sion. Based upon these we then assess the implications for the 
number, variability, and redshift evolution of the TeV blazars. 

In general there exist two main classes of AGNs that differ 
in their accretion mode and in the physical processes that dom- 
inate the emission. 

1. Thermal/disk-dominated AGNs: Infalling matter assem- 
bles in a thin disk and radiates thermal emission with a 
range of temperatures. The distributed black-body emis- 
sion is then Comptonized by a hot corona above the 
disk that produces power-law X-ray emission. Hence 
the emission is a measure of the accretion power of the 
central object. This class of objects are called QSOs or 
Seyfert galaxies and make up about 90% of AGNs. They 
preferentially emit in the optical or X-rays and do not 
show significant nuclear radio emission. None of these 
sources have so far been unambiguously detected by 
Fermi or imaging atmospheric Cherenkov telescopes be- 
cause the Comptonizing electron population is not highly 
relativistic and emits isotropic ally, i.e. there is no beam- 
ing effect that boosts the emission. 

2. Non-thermal/jet-dominated AGNs: The non-thermal 
emission from the radio to X-ray is synchrotron emis- 
sion in a magnetic field by highly energetic electrons that 
have been accelerated in a jet of material ejected from 
the nucleus at relativistic speed. The same population of 
electrons can also Compton up-scatter any seed photon 
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population either provided by the synchrotron emission 
itself or from some other external radiation field such as 
UV radiation from the accretion disk. Hence the SED of 
these objects shows two distinct peaks. The luminosity 
of these non-thermal emission components probes the jet 
power of these objects. Observationally, this leads to the 
class of radio-loud AGNs which can furthermore be sub- 
divided into blazars and non-aligned non-thermal dom- 
inated AGNs depending on the orientation of their jets 
with respect to the line of sight. 

There are no known sources above 600 GeV that correspond 
to AGNs with jets pointed at large angles (« 15° -40°, see 
Urry & Padovani 1995) with respect to the line of sight (for an 
example of a non-aligned AGN, NGC1275, that shows a very 
steep high-energy spectrum, emitting a negligible number of 
VHEGRs, see Mai'iotti & MAGIC Collaboration 2010). Hence 
we turn our attention to blazars, which can be powerful TeV 
sources. 

Blazars can further be subdivided into two main subclasses 
depending upon their optical spectral properties: flat spectrum 
radio quasars (FSRQ) and BL Lacs. FSRQs, defined by broad 
optical emission lines, have SEDs that peak at energies below 
1 eV, implying a maximum particle energy within the jet and 
limiting the inverse-Compton scattered photons mostly to the 
soft gamma-ray band. It is presumably for this reason that 
no continuous TeV component has been detected in an FSRQ 
(note, however, that TeV flares from FSRQs have been detected 
in two cases (MAGIC Collaboration et al. 2008b; Mose Mari- 
otti 2010)). 

In contrast, BL Lacs or Blazars of the BL Lac type (Mas- 
saro et al. 2009) can be copious TeV emitters. These are very 
compact radio sources and have a broadband SED similar to 
that of strong lined blazars, though lack the broad emission 
lines that define those. Depending upon the peak energy in the 
synchrotron spectrum, which approximately reflects the maxi- 
mum particle energy within the jet, they are classified as low-, 
intermediate-, or high-energy peaked BL Lacs, respectively 
called LBL, IBL, and HBL (Padovani & Giommi 1995; Abdo 
et al. 2010d). While LBLs peak in the far-IR or IR band, they 
exhibit a flat or inverted X-ray spectrum due to the dominance 
of the inverse-Compton component (see Fig 15 of Abdo et al. 
2010d, for a visualization of the SED of BL Lacs). The syn- 
chrotron component of IBLs peaks in the optical which moves 
their inverse-Compton peak into the gamma-ray band of Fermi. 
HBLs are much more powerful particle accelerators, with the 
synchrotron peak reaching into the UV or, in some cases, the 
soft X-ray bands. The inverse-Compton peak can then reach 
TeV energies (Ghisellini & Tavecchio 2008; Tavecchio & Ghis- 
ellini 2008; Abdo et al. 2010d).i^ 

In the gamma-ray band, the subclass of IBLs that emit VHE- 
GRs are almost indistinguishable from the HBLs, suggesting 
that the location of the synchrotron peak does not uniquely 
characterize the VHEGR emission from these sources (e.g., 
due to variations among individual blazars in the magnetic field 
strength within the synchrotron emitting region and the origins 

The source classes of HSP/ISP/LSP used in recent Fermi publications are 
very similar' to the commonly used HBL/IBL/LBL classes. Hence we identify 
these with each other, respectively, though minor differences may be found in 
the literature. Nevertheless, where we refer to the number counts observed by 
Fermi specifically, we will refer to the classes HSP/ISP/LSP in keeping with 
their notation. 

For the highest energies scattering occurs in the Klein-Nishina regime, 
which results in steeper inverse-Compton spectra. 



and properties of the seed photons that are ultimately Comp- 
tonized). Hence we identify HBLs and VHEGR-emitting IBLs 
with the single source class of TeV blazars. We note that there 
is presently no evidence for the hypothetical class of ultra- 
HBLs that were proposed to have a very energetic synchrotron 
component extending to 7-rays (Ghisellini 1999). If such a 
population of bright and numerous sources exists, Fermi should 
have seen it (Abdo et al. 2010c). The ultra-HBLs may have 
escaped detection from Fermi thus far by being either intrinsi- 
cally dim 7-ray sources or very rare objects (Costamante et al. 
2007). 

TeV blazars have a redshift distribution that is peaked at low 
redshifts extending only up to z = 0.7. This is most likely en- 
tirely a flux selection effect; TeV blazars are intrinsically less 
luminous than LBLs and FSRQs, with an observed isotropic- 
equivalent luminosity range of 10'*'* -2 x lO"**" erg s"', with the 
highest redshift TeV blazars also being among the most lumi- 
nous objects (see Figures 23 and 24 in Abdo et al. 2010c). That 
TeV blazars should be intrinsically less luminous than FSRQs 
is not entirely unexpected, however Ghisellini et al. (2009) 
have argued that the physical distinction between FSRQs and 
TeV blazars has its origin in the the different accretion regimes 
of the two classes of objects. Using the gamma-ray luminosity 
as a proxy for the bolometric luminosity, the boundary between 
the two subclasses of blazars can be associated with the accre- 
tion rate threshold (nearly 1 % of the Eddington rate) separating 
optically thick accretion disks with nearly Eddington accretion 
rates from radiatively inefficient accretion flows. The spectral 
separation in hard (BL Lacs) and soft (FSRQs) objects then 
results from the different radiative cooling suffered by the rel- 
ativistic electrons in jets propagating into different surround- 
ing media (Ghisellini et al. 2009). Hence in this model, TeV 
blazars cannot reach higher luminosities than approximately 
2 X lO'^^ergs"' since they are limited by the nature of ineffi- 
cient accretion flows that power these jets and by the maximum 
black hole mass, ^ 10'"Mo. 

5. 1.2. An Empirical Local (f>B( z,L) 

While we have attempted to place the observed TeV blazars, 
which we have identified with the HBLs and VHEGR-emitting 
IBLs, into the broader context of AGNs using the unified 
model, due to the substantial distinctions in accretion rate, 
emission properties, object morphology and geometry, it is not 
obvious that any of the properties of TeV blazars should be sim- 
ilar to those of AGNs more generally. Nevertheless, evidence 
for a simple connection between the two populations can be 
found in the similarity between their the luminosity functions 
(a fact we will exploit later in estimating the redshift evolu- 
tion of the TeV blazars). Here we define the luminosity for the 
purposes of defining (pB to be the isotropic-equivalent value as- 
sociated with emission between 100 GeV and 10 TeV. While 
this may be considered to be a VHEGR luminosity, because 
most TeV blazars are peaked within this band, this corresponds 
to the majority of the emission from these sources. 

The objects listed in Table 1 were chosen because they have 
well defined SEDs, based upon a combination of VERITAS, 
H.E.S.S., and MAGIC observations. These 28 sources have 
VHEGR spectra that are well fit by the form. 



dN / E 



(26) 



where /o is the normalization in units of cm' -^s-'TeV"'. The 
gamma-ray energy flux is trivially related to dN/dE by Fe = 
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Figure 5. Comparison between the luminosity-weighted quasar and TeV- 
blazar luminosity functions (L(f>Q{z,L) and L(I>b(z,L), respectively). The solid 
lines show the absolute L4>q (in comoving Mpc), while the dashed lines show 
L4iQ rescaled in magnitude by 2.1 X 10"^ and shifted to lower luminosities by 
a factor of 0.55. In both cases the black, orange, and red lines coiTespond to 
0g(O.l,L), (/)g(0.5,L), and 0g(l,L), respectively. The points and upper-limits 
show <^B of the HBL and IBL sources listed in Table 1 , assuming the relevant 
limits for PKS 1553-Hl 13 and PKS 1424-H240. Presented in the inset is the TeV 
source luminosity distance as a function of source luminosity for all of the 
sources in Table 1 with redshift estimates (including limits). The dotted line 
shows the distance-dependence of the flux limit we employ in the completeness 
connection. 

EdN /dE (X £■'"", from which we obtain a VHEGR flux, 



F = Eofo 



(•lOTeV 




/ dE 


/lOOGeV 





1-Q 



(27) 



and for sources with a measured redshift a corresponding 
isotropic-equivalent luminosity, L = AirDj^F, where Dl is the 
luminosity distance.'^ The resulting /o, Eq, a, F, and L are 
collected in Table 1. In addition we list the redshift, inferred 
comoving distance, and absorption-corrected intrinsic spectral 
index, defined at £0, obtained via 



d\nE 



a-TE[Eo{l+z),z]. (28) 



£0 



For high-redshift sources a can be substantially less than 2, 
implying that an intrinsic spectral upper-cutoff must exist. 

To produce 4>b, we must account for a variety of selection ef- 
fects inherent in the sample listed in Table 1 . The objects in Ta- 
ble 1 were originally selected for study for a variety of source- 
specific reasons, e.g., existing well known sources, extremely 
high X-ray to radio flux ratio in the Sedentary High energy 
peaked BL Lac catalog, hard spectrum sources in the Fermi 
point source catalog, and flagged as promising by the Fermi- 
LAT collaboration. In addition, the source selection suffered 
from the usual problems associated with surveys (e.g., schedul- 
ing conflicts with other targets, moon, bad weather, etc.). As a 

" Since the VHEGR spectra of TeV blazars typically are peaked above 
lOOGeV, this overestimates the luminosity by a factor of order unity. 



consequence, this sample is somewhat inhomogeneous. Nev- 
ertheless, in lieu of a less-biased sample, we will treat it as ho- 
mogeneous and correct for the selection effects were possible, 
focusing upon those due to the sky coverage and duty cycle of 
TeV observations, and those due to sensitivity limits of current 
imaging atmospheric Cerenkov telescopes. 

To estimate the sky completeness and duty cycle of this 
set of objects we rely upon the all-sky GeV gamma-ray ob- 
servations of HBL and IBL sources by Fermi (Abdo et al. 
2010c). Outside of the Galactic plane, Fermi observes 118 
high-synchrotron peaked (HSP) blazars and a total of 46 
intermediate-synchrotron peaked (ISP) blazars. Roughly half 
of the latter are likely to emit VHEGRs as indicated by their 
flat spectral index between 0.1 and 100 GeV (a < 2; see the 
spectral index distribution of Figure 14 in Abdo et al. (2009)). 
Of these potential 141 TeV blazars, only 22 have also been co- 
incidentally identified as TeV sources, whereas there are a total 
of33 known TeV blazars (29 HBL, 4 IBL). If these 141 sources 
are all TeV emitters, but have not been detected due to incom- 
plete sky coverage of current TeV instruments, then the selec- 
tion factor is rysei = 141 /33 = 4.3. In addition, the duty cycle 
of coincident GeV and TeV emission is T^duty = 33/22 = 1.5.'^ 
Finally, by excluding the Galactic plane for galactic latitudes 
h < 10°, this is an underestimate by roughly rj^ky = 1.17. 

We make a rough attempt to correct for the sensitivity limit of 
the TeV instruments, inferring a flux limit by fitting the upper- 
envelope of the TeV-blazar flux-luminosity distance distribu- 
tion. Somewhat surprisingly, despite the heterogeneous nature 
of the TeV observations, are remarkably well described by a 
single flux limit: 4.19 x 10~'^ergcm~^ s"'. This process and 
the associated limit are shown explicitly in the inset of Fig- 
ure 5. From this we obtain a maximum redshift, and therefore 
comoving volume, associated with each luminosity. We then 
construct (f>B by counting the number of TeV blazars per unit 
logjgL in each logarithmic luminosity bin and dividing by the 
comoving volume. The resulting (j>B, weighted by luminos- 
ity (and therefore showing the luminosity density in comoving 
units), is shown in Figure 5. It peaks at ~ 4 x 10'*'* ergs"', im- 
plying that, as expected, these objects are systematically dim- 
mer than most other AGN, and exhibits the broken-power law 
shape typical of AGN luminosity functions. More importantly, 
despite a handful of sources with redshifts ^ 0.5, the objects 
in Table 1 are all nearby, and therefore our (f>B corresponds to 
that in the local Universe, i.e., z ^ 0.1. Based upon this, the 
inferred present-day TeV-blazar luminosity density is roughly 
1 X lO^'^ergs'^Mpc'l 

Also shown in Figure 5 is the L(j)Q obtained by Hopkins et al. 
(2007). After rescaling LipQ to lower luminosities (0.55) and 
lower luminosity densities (2.1 x 10""'), it provides a remark- 
ably good fit to our L(j>B (reduced-x^ of 0. 13 with 3 degrees of 
freedom), i.e., we find 



/)b(0.1,L)~3.8x 10">e(0.1,1.8L), 



(29) 



where an explicit expression for (f>Q is given in Appendix B. 
While there is considerable uncertainty in Lips, especially at 
low luminosities, it clearly does not fit L(f>Q at higher z- This 
suggests two immediate conclusions: 

'** Note that we implicitly assume that the luminosity distribution of ob- 
served TeV sources reflects the true distribution after accounting for flux in- 
completeness, and thus constant correction factors (independent of luminosity) 
can be used to estimate the true distribution. We show that this assumption is 
fulfilled in Figure 5. 
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Figure 6. The normalized number of TeV blazai's (presumably the hai'd 
gamma-ray blazars comprised of hard ISPs and all HSPs) that are expected 
to have been observed by Fermi as a function of redshift. The dashed lines and 
solid histograms show the ideal and binned distributions for a = 1.95, though 
neither varies significantly with a. The total number of objects is sensitive to 
the spectral index, and shown as a function of a in the inset. Different col- 
ors con'espond to different normalizations between the lOOGeV-lOTeV and 
lOOMeV-lOOGeV luminosities: r^mm = 0.78 (red), 1.6 (green), and 3.1 (blue), 
i.e., the former is the latter luminosity multiplied by »)n,iii)- For reference, the 
normalized number of hard gamma-ray blazars (HSPs and half of the ISPs) 
observed by Fermi are shown by the black circles, and the black square in the 
inset shows the average spectral index and total number observed, with hori- 
zontal error bars giving the l-cr range. 

1. The bolometric output of TeV blazars and quasars are 
regulated by similar mechanisms, presumably accretion, 
despite the large difference in luminosity and the details 
of the emission processes between the two populations. 

2. TeV blazars and quasars are contemporaneous elements 
in a single AGN distribution; specifically, TeV-blazar ac- 
tivity does not lag that of quasars. 

5.1.3. Extending 4>b(z,L) to high z 

Based upon the strong similarities between (pB and (pQ, and 
the associated implications, we make the conservative theoret- 
ical assumption that the redshift evolution of the TeV blazars 
follows that of quasars. That is, we suppose that Equation (29) 
holds at all z- This implies that the integrated comoving TeV- 
blazar isotropic-equivalent luminosity density is given by 

/oo 
L0B(z,L)t/logioL = 77BAe(z), (30) 
oo 

where the constant of proportionality, ry^ ~ 2.1 x 10"^, is then 
set by the comparison between (pB and 0q at z = 0.1. As a con- 
sequence, in comoving units, the TeV-blazar luminosity den- 
sity would be roughly an order of magnitude larger at z ~ 2, 
~ 1 X lO^^'^ergs-'Mpc-^ 

5.2. Fermi TeV Blazar Counts and the Local Evolution of the 
Blazar Luminosity Function 

The redshift distribution of the Fermi BL Lac sample (i.e., 
the First LAT Catalog, ILAC; Abdo et al. 2010c) is peaked at 
z < 0.2 and falls rapidly thereafter Inasmuch as the Fermi HSP 
and ISP counts directly probes the low-z (pB, it may appear that 
they are inconsistent with the rapidly evolving (pB described 
in the previous section. Indeed, an analysis of the first three 



months of Fermi observations, suggested that the BL Lac pop- 
ulation does not grow substantially with redshift (Abdo et al. 
2009). However, as seen in Figure 6, this is not necessarily the 
case. 

Assuming, as we have, that the number of TeV blazars is 
equal to the the number of hard Fermi gamma-ray blazars, cor- 
responding to ISPs with a < 2 (comprising roughly half of that 
class) and HSPs, the expected number observed inside of a 
given redshift is: 

Mb{z}= rdz'^4nDl 
Jo dz' 

X / °'° °" (l+z')VBfe',i)dlogioi, (31) 

where D = j cdt' = J cdz' / [//(z')(l +z')] and Da = Dl/{1 + 
z)^ are the proper and angular diameter distances, respectively, 
Ln,ax — 2 X 10'*^ ergs"' is the intrinsic upper-cutoff of the TeV 
blazar isotropic-equivalent luminosity function, and 

imm(z)=^min47rD2(z)(l+z)""' 



~4x 10^" rj, 





'DdzY 







is lower-limit set by the Fermi flux limit (see Figure 23 
of Abdo et al. 2010c, and surrounding discussion). The 
factor 77,nin is a correction relating the lOOGeV-lOTeV 
isotropic-equivalent luminosities we employ to define (pB to 
the lOOMeV-lOOGeV luminosities used by Fermi to define 
the flux limit. The dependence upon a arises from the lim- 
ited spectral coverage of both definitions, though here we fix 
a = 3 for all objects based upon the sources that dominate the 
TeV flux at Earth. Note that since the (pQ from Hopkins et al. 
(2007) diverges at small L, the lower-luminosity cutoff is crit- 
ical to getting both the total number of TeV blazars and the 
shape of their redshift distribution correct. The resulting Mb is 
shown in Figure 6 for three different choices of t^nun'- 0.78, 1 .6, 
and 3.1, of which Tymin = 1.6 is most similar to the ILAC hard 
gamma-ray blazars. 

Generally, our (ps does an excellent job of reproducing the 
overall number of ILAC hard gamma-ray blazars and the dom- 
inance of nearby objects in their redshift distribution. This is 
despite the strong redshift evolution of (pB implied by its rela- 
tionship with (pQ. The reason for this is the flat distribution at 
low L, the steep drop-off at high L (a result of which is that the 
shape is only marginally sensitive to the cutoff at L^ax) and the 
rapidly growing L,nin due to the fixed flux limit. 

However, we note that the comparison between (pB and the 
ILAC hard gamma-ray blazar statistics assumes that those 
with measured redshifts, comprising roughly half of the ILAC 
hard gamma-ray blazar sample, are representative of the hard 
gamma-ray blazar population as a whole. This appears not 
to be the case; as we discuss in detail in Appendix C, it is 
clear that the ILAC HSPs with and without measured red- 
shifts are not drawn from the same underlying a-distribution. 
Based upon a similar analysis for ILAC BL Lacs generally, 
Abdo et al. (2010c) argued that the objects without redshifts 
are more consistent with z > 0.5 population. However, this 
conclusion does not extend to HSPs, for which there are only 
three sources in the clean ILAC HSP sample with z > 0.5. 
Rather, in Appendix C, we show that the ILAC HSPs without 
redshifts are distributed both in spectral index and flux much 
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more similarly with nearby HSPs (z < 0.25) than more distant 
HSPs (z > 0.25). If this is because these sources are intrinsi- 
cally under-luminous, nearby objects, it could marginally im- 
prove the already remarkable comparison between the number 
of objects implied by our <f>B and those observed. However, be- 
cause the number of predicted ILAC hard gamma-ray blazars 
is strongly dependent upon the flux limit, even a marginal in- 
crease in either 0^ at low luminosities or the effective flux limit 
results in a dAfg /dz that is more strongly weighted at low z, po- 
tentially allowing even more dramatically evolving luminosity 
functions. For this reason we conclude that our (f>B is broadly 
consistent with the ILAC hard gamma-ray blazar distribution, 
and that Fermi does not, at present, exclude a rapidly evolving 
<f)B in the recent past. 

5.3. The TeV Blazar Contribution to the Fermi EGRB 

While the statistics of the ILAC hard gamma-ray blazar sam- 
ple probes the evolution of at z < 0.5, the Fermi EGRB 
is sensitive to TeV Blazars at high z. The contribution to the 
EGRB from TeV blazars given our (pB is shown in Figure 7. Be- 
cause ICCs are suppressed, computing the EGRB flux requires 
only summing the individual intrinsic GeV spectra of the TeV 
blazars. For simplicity, to estimate the TeV blazar contribution 
to the EGRB we assume all TeV blazars have identical intrinsic 
spectra, given by a broken power law. 



Ff = LFf oc 



1 



'E/E,Y^-\(E/E„Y-' 



(33) 



/• lOTeV 

where the normalization is set such that the j^QQ^^yFEdE = 1 

(i.e., Jjylj^^ Lf^iiii = L, corresponding to the luminosity used 
to define (pB, and thus determine ?/b), Eb is the energy of the 
spectral break, and a/, < a are the low and high-energy spectral 
indexes, respectively.'^ With this, after performing the integral 
over L, the EGRB flux is then 



dE^ ' ^ 47ri; dz' AttDI 



Hi 

47r 



= , cAq(z') 
dz „ . 



(34) 



H(z'){l+zT 



where E' = E(l+z'), and Aq = Aq(1 is the physical quasar 
luminosity density. We choose fiducial values for the spectral 
parameters of q;l = 1.33, isi = 1 TeV, and a = 3, typical of the 
TeV blazars (see Table 1 and Ghisellini 2011). Finally, since 
Fermi resolves individual gamma-ray blazars with isotropic- 
equivalent luminosities ~ 10"*^ ergs"', roughly the location of 
the peak in 4>b, for z < 0.25, we construct the EGRB from 
sources at larger redshift. 

If Eh is sufficiently high (> 2 x 10" GeV) the EGRB is re- 
markably insensitive to the particular values of Eh and a (see 
the bottom two panels of Figure 7). This is a direct result of 
the annihilation of the VHEGRs upon the EBL, effectively re- 
moving the relevant portion of the blazar SED from the EGRB. 
More important is the low-energy spectral index. However, for 
values of at that are consistent with the Comptonization mod- 
els for the VHEGR emission, the anticipated EGRB is consis- 

" Note that 4>b(7.,L) was detemiined assuming an unbroken spectrum, and 
thus in this case would over estimate the total flux from the TeV blazars. Nev- 
ertheless, here this is not accounted for, i.e., we renormaUze the broken power- 
law spectrum to our previous power-law estimate of the TeV blazar luminosity 
density. As a consequence, our estimates of the EGRB are also over estimates. 
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Figure?. Contribution of TeV blazars to the EGRB. The solid lines show 
the contribution from TeV blazars beyond z, = 0.25, roughly the set for which 
Fermi can no-longer resolve them. The dashed Hues show the contribution 
from all TeV blazars, including those which Fermi should have been able to 
detect individually and to remove. The dotted Hues show the contribution from 
all TeV blazars when the annihilation upon the EBL is neglected. Different 
curves correspond to different parameters of the intrinsic spectrum assumed for 
the blazar populations. Top: Variations in the low-energy spectral index (top to 
bottom, ql = 1.83, 1.67, 1.5). Middle: Vaiiations in the break energy (top to 
bottom, E], = 500 GeV, 1 TeV, 2 TeV). Bottom: Variations in the high-energy 
spectral index (top to bottom, a = 3.5, 3, 3.5). Unless otherwise specified, the 
remaining parameters are = 1.63, Ei, = 1 TeV, and a = 3. Finally, the Fermi 
measurement of the EGRB reported in Abdo et al. (2010a) is shown by the 
blue points. Note that in our scenario, the integral of the difference between 
the unabsorbed (dotted) and absorbed (dashed) energy fluxes, which dominates 
the total energy budget of these sources, has been dissipated into heat within 
the IGM. 
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tent with the Fermi result. Thus, we find that despite our dra- 
matically evolving 0^, we are are able to satisfy the limits im- 
posed by the Fermi EGRB. Moreover, for reasonable spectral- 
parameter values it is possible to accurately reproduce both the 
magnitude and shape of the high-energy EGRB (i.e., at ener- 
gies > lOGeV).^*^ At first this may appear in conflict with other 
studies that have performed more sophisticated analyses of the 
gamma-ray emission from blazars (e.g., Narumoto & Totani 
2006; Kneiske & Mannheim 2008; Inoue & Totani 2009; Ven- 
ters 2010) that have typically found that such a strongly evolv- 
ing 4>B over-produces the Fermi EGRB. However, this is not 
the case for three reasons. 

First, the AGN populations typically considered in computa- 
tions of the EGRB include the FSRQs and LBLs, which domi- 
nate at low energies. The TeV blazars of interest here are only 
the high-energy tail of the blazar distribution, and it is these 
alone that we fix to the quasar luminosity function (indeed, it 
is for only these objects that the arguments surrounding Figure 
5 have been made). Moreover, the TeV blazars themselves are 
generally dimmer, both bolometrically and within the Fermi 
LAT energy band (200MeV-100GeV), than the LBLs and FS- 
RQs, which have typical isotropic-equivalent luminosities of 
3 X 10"*^ ergs"' (Narumoto & Totani 2006). 

Second, and perhaps most importantly, the suppression of 
the ICCs means that the VHEGR emission is lost to heat in the 
IGM, not reprocessed to below lOOGeV. Since the TeV blazar 
flux is presumed to be dominated by emission above lOOGeV, 
the non-radiative dissipation of this emission substantially re- 
duces the impact of the TeV blazars upon the EGRB. For all of 
the spectral parameters we considered the un-absorbed fluxes 
exceed the 30GeV EGRB by a significant margin, implying 
that the suppression of the ICCs is crucial to bringing the TeV 
blazar contributions in line with the Fermi EGRB. Thus, the 
lack of the ICCs generally appears necessary to reconcile the 
blazar and quasar luminosity functions. 

Third, the VHEGR spectra from the most luminous TeV 
blazars necessarily peaks near or above TeV energies. For a 
number of the sources in Table 1, the intrinsic VHEGR spec- 
trum (adjusted for annihilation on the EBL) is inverted, indicat- 
ing that this break is above a TeV. The brightest TeV blazar in 
Table 1, Mkn 421, has an inverted spectrum below lOOGeV 
(Abdo et al. 2010c), implying Eh ^ 500 GeV. The second 
brightest, lES 1959H-650, is consistent with being flat below 
100 GeV, and appears to peak near ITeV when in the high- 
state (Daniel et al. 2005). Similarly, lES 2344H-5 14 peaks near 
~ 300 GeV (Albert et al. 2007c) and Mkn 501 is inverted be- 
low 100 GeV, implying a peak above that value (Abdo et al. 
2010c). Thus the sources likely to dominate the VHEGR 
background, and by extension the EGRB, all turn over near 
lOOGeV-1 TeV. This is expected if the the VHEGR emission 
arises from inverse-Compton scattering the synchrotron bump 
(see, e.g., Ghisellini 2011, and references therein). The con- 
sequence of the spectral break is the suppression of the contri- 
bution from the TeV blazars to the EGRB below Eb- Were the 
TeV blazar contribution to the EGRB dominated by sources 
with Eb <2x 10^ GeV it would exceed the Fenni EGRB at 
E < lOGeV. 

A more complete analysis of the EGRB should include a 
variety of spectra, including a distribution of break energies, 
smoothly connecting the HBL, IBL and LBL populations. 

Note that since we are only presenting the EGRB from the TeV blazars, 
we necessarily neglect the contributions from the FSRQs and LBLs which 
dominate the EGRB at low energies. 



However, this is beyond the scope of this paper, which is pri- 
marily concerned with the fate of the VHEGR emission, absent 
in the LBLs and FSRQs. Nevertheless a recent comprehen- 
sive model, which includes the contributions of the FSRQs and 
BL Lacs observed by Fermi, as well as starburst galaxies, has 
had considerable success fitting the Fermi EGRB with an even 
more extremely evolving (f>B', fixing it to the luminosity func- 
tion of radio galaxies (Cavadini et al. 2011; Stecker & Ven- 
ters 2011). Of particular relevance here is that Cavadini et al. 
(201 1) and Stecker & Venters (201 1) explicitly ignored the po- 
tential contributions of ICCs. Unlike their analysis, however, 
we find that the TeV blazars are capable of reaching the high- 
est Fermi bands despite the annihilation upon the EBL. 

6. CONCLUSIONS 

The cold, highly anisotropic beams of ultra-relativistic e^ 
pairs produced by the annihilation of VHEGRs upon the EBL 
are unstable to plasma beam instabilities. More importantly, 
for a wide range of parameters relevant for the observed TeV 
blazars these instabilities may be capable of isotropizing, and 
potentially extracting the kinetic energy of, the pairs at a rate 
orders of magnitude faster than inverse-Compton scattering. 

This has far reaching consequences for efforts to constrain 
the IGMF using empirical limits upon the GeV emission from 
known TeV sources. Typically, ~ 300 yr after the onset of 
TeV emission the pair beam density has grown sufficiently for 
plasma beam instabilities to dominate its evolution, randomize 
the beam, and potentially suppress the inverse-Compton sig- 
nal upon which the IGMF limits are based. Note that due to 
the beam disruption by the instabilities, this occurs even if the 
plasma instabilities do not ultimately cool the pairs. As a con- 
sequence, the present constraints upon the IGMF, obtained by 
the non-observation of an inverse-Compton GeV bump in the 
spectra of bright TeV blazars are inherently unreliable. Nev- 
ertheless, the sudden appearance of a TeV-bright blazar or in- 
trinsically transient sources (e.g., GRBs) provide a means to 
temporarily avoid the consequences of plasma beam instabili- 
ties during the growth of the pair beam. Alternatively, observ- 
ing particularly dim sources, L < lO'^^ergs"', limits the beam 
density directly, again avoiding the complications imposed by 
plasma processes. Finally, the presence of these plasma insta- 
bilities in the pair beams of TeV blazars, which manifest them- 
selves through their impact on the IGM (see Paper II, Paper III, 
and Puchwein et al. (2011)), implies the most stringent upper 
limit to date on the IGMF: B < W'^G. 

If the plasma instabilities can efficiently convert the pair 
beam kinetic energy into heat in the IGM, as we anticipate 
based upon existing numerical simulations and the arguments 
in Section 3.5, they would necessarily suppress the develop- 
ment of ICCs, and thus prevent the reprocessing of the VHEGR 
emission from bright sources to GeV energies. The lack of 
ICCs, independent of the mechanism that facilitates the local 
dissipation of the pair kinetic energy, would greatly weaken the 
constraints upon the evolution of the blazar population derived 
from the unresolved EGRB measured by Fermi. By introduc- 
ing a spectral break near 1 TeV and eliminating the reprocessed 
VHEGR emission, we find that the Fermi EGRB is consistent 
with a TeV blazar (and therefore presumably hard gamma-ray 
blazar) luminosity function fixed to that of quasars, normalized 
by comparing objects in the local Universe (z ~ 0.1), and mo- 
tivated by the remarkable similarity between them in the local 
Universe. This conclusion is relatively insensitive to the par- 
ticular parameters governing the VHEGR spectra of the TeV 
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blazars, requiring only that the VHEGR emission is produced 
via inverse-Compton scattering. For a wide range of spectral 
parameters, we are able to match the magnitude and shape of 
the Fermi EGRB at high energies, with the low-energy compo- 
nent presumably arising from the FSRQs and LBLs. This (t>B, 
and perhaps an even more rapidly evolving luminosity func- 
tion, is also consistent with the observed redshift distribution 
of the ILAC HSP and hard-ISP sample. 

Matching the high-energy EGRB (above lOGeV) requires 
more TeV blazars than are currently observed, though a 
comparable number to those inferred once the sky-coverage 
and GeV duty-cycle completeness corrections are included. 
Based upon these factors, and our rough estimate of the 
EGRB, we predict that upcoming surveys performed with next- 
generation Cerenkov telescope arrays (see, e.g., CTA Consor- 
tium 2010) should find roughly 2 x 10~ sources above 4.2 x 
10~'~ erg s"' cm"- (our estimate of the effective flux limit of cur- 
rent imaging Cerenkov telescopes), and a handful of additional 
sources comparable to the brightest TeV blazars observed. In 
addition, based upon the current number of known TeV blazars 
and our estimate of 4)b{z,L), the improved anticipated sensi- 
tivities of these instruments, 5-10 time larger than current ar- 
rays, should result in the detection of 1.5 x 10^-3 x 10^ addi- 
tional TeV blazars, with median luminosities ^3x10"*^ erg s"' . 
These should allow more precise estimates for their gamma- 
ray SEDs and a better characterization of <j>B{z,L), especially 
for low-luminosity objects. 

Unlike inverse-Compton cooling, the plasma beam instabil- 
ities deposit the energy locally, heating the IGM. Moreover, 
the homogeneity of the EBL, and the weak dependence of the 
plasma cooling rates upon the IGM density, result in a uniform 
volumetric heating, in clear contrast to either photoionization 
heating or mechanical feedback from AGN. While we shall 
defer a detailed discussion of the consequences of this heat- 
ing to Papers II, III, and Puchwein et al. (2011) here we note 
that this unusual heating prescription naturally explains a num- 
ber of heretofore outstanding questions, including the inverted 
equation-of-state (temperature-density relation) for low density 
regions in the IGM (Paper II), the suppression of dwarf galax- 
ies and their histories, the segregation of galaxy clusters and 
groups into cool core and non-cool core populations (Paper III), 
and the quantitative properties of the high-redshift Lya forest 
(Puchwein et al. 2011). As a consequence, despite the fact that 
our estimates of the plasma cooling rates are limited to the lin- 
ear regime (though with some numerical support), there are a 
variety of observational reasons to believe that plasma cooling, 
or an analogous mechanism, does in fact dominate the evolu- 
tion of the ultra-relativistic pair beams. 
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APPENDIX 
A. INSTABILITY GROWTH RATES 

Here we compute the growth rates for the various plasma 
instabilities discussed in the text within the context of relativis- 
tically moving pair plasmas. In all cases we make use of the 
kinetic theory description of the underlying plasmas. This nec- 
essarily assumes that the plasma, and in particular the beam 
plasma, is sufficiently dense that it is well described by a distri- 
bution function on the relevant scales. This is equivalent to re- 
quiring that many particles within a characteristic energy range 
be present on the plasma scale of the IGM, i.e., the conditions 
laid out in Section 3.1. In our analysis the Vlasov equation will 
play a central role, which owing to the relativistic nature of the 
calculation, we express in terms of the canonical pair, x and p, 
making use of the Lorentz invariance of the electron/positron 
distribution functions /T(x,p) and the phase-space volume el- 
ement cPxcPp. In what follows we set c = 1 unless otherwise 
specified. 

A.l. Relativistic Pair Two-Stream Instability 

The two-stream instability arises due to the excitation of 
negative-energy electrostatic waves in the beam and target plas- 
mas. These waves carry away both the energy and momentum 
of the beam. Specifically, we compute the growth rate of the 
electrostatic wave moving in the direction of the beam in the 
absence of a background magnetic field^'. In this case, the 
Vlasov equations for the electrons and positrons are; 



Dt 

or 



op 



dr 

+eE.^ = 0, 
Dt 9p 



(Al) 



where D/Df = d/dt+\- V is the Stokes derivative and E is the 
net electric field. Linearizing these and Fourier transforming in 
t and X gives 



-/(^-k.v)/r-eEi.^ = 
-i{uj -k-y) ft + eEi--^ = 0, 



(A2) 



where are the perturbations to the electron and positron 
distribution functions, Ei is the electrostatic wave field and 
we have assumed that the background distributions, f^^, are 
isotropic. As a result we have 



±ie 



-El 



dp 



(A3) 



At this point we may compute the dielectric tensor associ- 
ated with the plasma response, however it will suffice to con- 
sider Gauss's law. Thus we now compute the perturbed charge 

^' Where the beam energy density is lai'ge in comparison to that associated 
with the background magnetic field this is well justified. 



THE COSMOLOGICAL IMPACT OF BLAZAR TEV EMISSION I 



19 



density: 



Pi=y {eft-ef;)d'p 



1 r^/o.Muv 



w-k-v \ dp dp 

a 1 



-El - 



dp w-k- V 
(k-k vv) 



(A4) 



7(w-k-v) 



It is now necessary to specify the We idealize the tar- 
get ionic plasma as cold and the pair beam plasma as mono- 
energetic, yielding 



f^ = n,5\p)+'^5\p-p,) and /J = y <S'(p-Pi) , (A5) 

where n, and «i are the lepton number densities in the target 
and beam, respectively. After performing the trivial integrals 
we then obtain 



ie^^ ( k k-k-ViVi 
pi = — El • n, — +nh— ; "2 



— E| -k -ir + 



rib 



^2 ■yl(uj-kVh)- J ' 



(A6) 



were in the final equality we used the fact that k || v/,. From 
Gauss's law we have ik ■ Ei = A-irpi, and therefore 



1 



■0, 



(A7) 



where Wp, = Ane^yi,/me and tjJ^b^ A-ire^rib/me are the plasma 
frequencies associated with the target and beam plasmas. 

This explicitly provides the dispersion relation, quadratic in 
u! (one electrostatic wave traveling in each direction for each 
plasma). When rib = 0, cj = ujpj. When iib ^ n,, as is the case 
of interest here, we may solve the dispersion relation perturba- 
tively. We do this by setting lu = ujp,{l+ri), with r/ <C 1, which 
gives: 



27J- 



= 0, 



(A8) 



which is no longer independent of k. When \ujpj-kvb\ ^ 
|?7aj/)., |, rj is real and thus there is no instability. On the other 
hand, where k ~ up^t/vb, we have 



27?- 



P.b 



3 2 T 

Tb'^pjT 



^0, 



and therefore ry' = n/j/2«,7^. This has three solutions: 



r] = 



1 / Hb 



lb \2ni 



1/3 



1 V3 I V3 
1, 1, — -I- ; 

2 2' 22 



(A9) 



(AlO) 



the first of which is oscillatory, the second is decaying and the 
third is growing with timescale 



^(uj) = Q{u!p,,r]) 



%/3 / Hb 
2 \2n, 



1/3 



lb 



(All) 



This differs from that associated with non-relativistic, ionic 
beams only by the factor of 1 / 7/,, arising due to time-dilation 
within the beam. 

Note that since the energy within the electrostatic wave is 
proportional to E^, the rate at which energy is removed from 
the beam is Fxs = 23 (oj). 

A.2. Relativistic Pair Weibel Instability 

The relativistic version of the Weibel instability has been 
discussed in some detail in the literature for the case of 
anisotropic, though symmetric beams (Weibel 1959; Fried 
1959; Yoon & Davidson 1987; Medvedev & Loeb 1999). Here 
we consider a the case of a dilute pair beam incident upon a 
much denser target plasma. This is essentially identical to the 
two-stream instability discussed previously, though coupling 
instead to an electromagnetic mode with k _L v;,. 

Again we begin with the linearized Vlasov equations for the 
electrons and positrons: 



^-.(E:+vxBO-^=0 
Dt op 

-^ + e(Ei+vxBi).^=0, 
Dt op 



(A12) 



where we have Bi = k x Ei/w from Faraday's law. Fourier 
transforming in t and x and solving for /j^ gives 



^_^/e[Ei+vx(kxEi)/c^] df^ 



/i+ = ±- 



oi-k- V 

= ±!f fEi + J^kV^. 
Ul \ w-k-v / dp 

The associated perturbation to the current is 



dp 



(A13) 



h = j {eft-ef;)yd'p 



dp 



v-E, 



A. 



1 + 



kv-i-vk 



• V 



w - k • V (a; - k • v)2 

(A14) 

in which we've defined /o = /o +/o ■ 

Choosing the given by Equation (A5) and assuming k || 
\b II El we recover the standard two-stream instability. For our 
purposes here, however, the computations may be substantially 
simplified by boosting into a frame in which the target plasma 
is not at rest. In this frame we have /o given by, 

fo = n',S\p'-p',) + n'bS\p'-p'b), (A15) 

with p', p^ II where we will take care at the end to properly 
relate all of these quantities to their target-frame analogs. In 
particular, we choose this frame such that «,v,' = «/,v^, where n, 
and lib are the proper densities of the target and beam plasmas, 
respectively. If Vb is the beam plasma velocity in the target (or 
lab) frame, this gives 



l-VbV', 



(A16) 



where £_=nb/ fit and is in our case much less than unity. Note 
that this is only the center of momentum frame if ^ = 1 (for 
which v'l = v^, the case most commonly discussed). If we 
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choose k' _L this choice of frame causes the terms linear in 
in Equation (A14) to vanish identically, yielding 



Ji 





4) 


A7/ 


lb) 



1, 



{All) 

Using the inhomogeneous wave equation, obtained from Fara- 
day's and Ampere's law. 



and choosing Ej || produces the dispersion relation 



(A18) 



It 



"^p,b 
it 




■'Rh^'b 

ib 



= 0. (A19) 



At this point we make use of the fact that in our case ^ ^ 1, 
which allows a perturbative solution of Equation (A16); v', = 
^Vi + 0(C^) and therefore vj, = Vfc(l-0 + C(^^)- Furthermore, 
Lo' = ijo + C'(^^) and k' = k for the geometry under considera- 
tion. As consequence, a;^,vf /7,' ~ ujp^,£,^vl = ^cop.bvl/lh < 

^p.hvihb and <jJp^,h? - = C^^l,b > ^'p,bhl'b^ which im- 
plies that 



UJ 



""Rb^'b 



lb 



(1-20 = 0. 



This has a purely imaginary solution: 



k2 + 



4kW„,vl(\-20 



(k^ + ioD^ji, 



1/2 



■1. 



(A20) 



(A21) 



For k ^ u!pj this rises linearly with k, saturating for k > wpj ^ 
i^p,bvl/ y/^, giving the growth rate 



5M: 



i^P.bVb 



(A22) 



Given that for our application ^ <C 1 , we will drop terms first 
order in ^ as well. Note that the plasma frequency that en- 
ters into the growth rate is that of the beam, which is much 
lower than that associated with the target plasma. Neverthe- 
less, the scale of the rapidly growing perturbations is limited to 
that associated with the target plasma, which are general much 
smaller than that associated with the perturbations in the beam 
alone. Again, the rate at which energy is sapped from the beam 
is then Fw = 23 (w). 

B. AN EXPLICIT EXPRESSION FOR THE QUASAR LUMINOSITY 

FUNCTION 

In the interests of completeness, here we reproduce the 0q 
from Hopkins et al. (2007), corresponding to the "Full" case in 
that paper, that we employ. See Hopkins et al. (2007) for how 
this (f)Q was obtained, and caveats regarding its application. 

The form of (f)Q is assumed to be a broken power law: 



<Pq{z,L) = 



4>* 



(Bl) 



[L/L,(z)]'/'fe) + [L/L,(z)]T2fe) ' 

where the location of the break (L*(z)) and the power laws 
(71 (z) and 72 (z)) are functions of redshift. These are given by, 

log 1 L* (z) = (log 10 ) y + ^z,, 1 ^ + k^2^- + 
72(z) = 272,0 (l0^-'2 '« + 10*~'2 ^«)"' 



(B2) 
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Figure 8. Cumulative (top) and differential (bottom) number of clean Fermi 
HSPs as a function of measured spectral index. Different line colors corre- 
spond to different subpopulations of the Fermi HSP sample: all HSPs with 
measured redshifts (green), HSPs with measured redshifts > 0.25 (red), HSPs 
with measured redshifts < 0.25 (blue), and HSPs without redshift measure- 
ments (black). In the bottom panel, error bai's denote the Poisson uncertainty 
only. 



where 



e = log 



1+z 



(B3) 



The values of the relevant parameters are given in Table 2, and 
where any estimate of the uncertainty is made we assume these 
are independently, normally distributed. 

C. INFERRING THE REDSHIFT DISTRIBUTION OF THE FERMI 
HSPS WITHOUT MEASURED REDSHIFTS 

Of the 118 HSPs detected by Fermi and reported in the First 
LAT AGN Catalog (ILAC, Abdo et al. 2010c, , there called 
HSPs), 113 are members of the clean sample (meaning that 
there are no ambiguities surrounding their detection), and of 
these only 65 have measured redshifts. This is a somewhat 
larger fraction that for ILAC BL Lacs generally, though this 
leaves nearly half of the HSP population with their redshifts 
undetermined. Based upon comparisons between the spectral 
index and flux distributions between the HSPs without redshifts 
and those at various redshift ranges, here we argue that these 
are likely to be located nearby. 

Already, based upon comparisons between the spectral in- 
dex distributions (SIDs) of the ILAC BL Lacs at large, it is 
clear that the objects with and without measured redshifts are 
not drawn from the same underlying population (Abdo et al. 
2010c). Similarities between the SIDs of the unknown-z ob- 
jects and the z > 0.5 subset of those with measured redshifts, 
Abdo et al. (2010c) has suggested that the BL Lacs without 
redshifts may be biased towards higher redshifts. However, we 
note that there are only three HSPs with z > 0.5, and thus this 
conclusion is relevant for ISPs and LSPs only. 

The SIDs of the unknown-z HSPs is shown in Figure 8, com- 
pared with the SIDs of HSPs with redshifts (cf. Figure 22 of 
Abdo et al. 2010c). The Kolmogorov-Smirnov (KS) proba- 
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Table 2 

Parameters of the Quasar Luminosity Function from Hopkins et al. (2007) 



Normalization 



71 



72 



log|( 



-4.825 ±0.060 



(logioi*)o 

kL.3 



13.036 ±0.043 
0.632 ±0.077 
-11.76±0.38 
-14.25 ±0.80 



71.0 



0.417 ±0.055 
-0.623 ±0.132 



72.0 
^72,2 



2. 174 ±0.055 
1.460 ±0.096 
-0.793 ±0.057 



s orcomoving Mpc 
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Figure 9. Cumulative (top) and differential (bottom) number of clean Fermi 
HSPs as a function of integrated flux between IGeV and lOOGeV. Different 
line colors correspond to different subpopulations of the Fermi HSP sample: all 
HSPs with measured redshifts (green), HSPs with measured redshifts > 0.25 
(red), HSPs with measured redshifts < 0.25 (blue), and HSPs without redshift 
measurements (black). In the bottom panel, error bars denote the Poisson un- 
certainty only. 



Figure 10. Redshift and spectral index-flux distribution of the clean Fermi 
HSP sample. HSPs with measured redshifts are shown by the red circles, with 
the circle size hnearly related to ; (ranging between z = and 0.7, i.e., large z 
is denoted by large points and small z by small points). HSPs without redshifts 
are shown by the blue triangles. For reference, the dotted black hues show 
hues of constant redshift for a source with luminosity 2 X lO"** ergs"' between 
lOOMeV and lOOGeV (the definition of L.^ in Abdo et al. 2010c, note the 
difference with the definition of F35). 



Table 3 

Comparison Between HSPs With and Without Measured Redshifts 



Redshift Spectral Index Flux between 1 GeV- 100 GeV 

Population Pks ("> ^'ks (logio (^35/cm"-s"')> 

All HSPs with z's 0.04 2.84 ±0.03 0.07 -8.89 ±0.01 

HSPs with z > 0.25 0.01 2.83 ±0.05 0.007 -9.06 ±0.03 

HSPs with z < 0.25 0.45 2.84 ±0.03 0.42 -8.78 ±0.02 

HSPs without z's — 2.91 ±0.03 — -8.86 ±0.02 



bility that these are drawn from the same parent population is 
0.04, indicating that this is unlikely at the 2-a level. In addi- 
tion, we show the SIDs of HSPs with z < 0.25 (blue line) and 
z > 0.25 (red line), with corresponding KS probabilities of 0.45 
and 0.01, respectively (these are collected in Table 3). Thus, in 
contrast to the ILAC BL Lac sample at large, the HSPs with- 
out redshifts appear to have SIDs that are strongly inconsistent 
with the population observed to have high redshifts, and in- 
distinguishable from those at low redshifts. Nevertheless, the 
HSPs with unknown-zs still tend to be softer than those with 



measured redshifts in either range. 

A similar analysis may be performed upon the reported flux 
measure, ^35, corresponding to the flux between IGeV and 
100 GeV. Figure 9 shows the flux distributions (FDs) of the 
unknown-z, all measured z, z > 0.25, and z < 0.25 HSP popula- 
tions, and is analogous to Figure 8. As before, the HSPs with- 
out redshifts have FDs that are marginally inconsistent with 
being drawn from the same parent population as the complete 
set of HSPs with measured redshifts, having a KS probability 
of 0.07. The inconsistency with the high-z HSP FD is even 
more striking than for the SIDs, with a KS probability less than 
0.007, i.e., that the high-z FD of the unknown-z and z > 0.25 
HSPs are from the same distribution is excluded. However, 
again, we find that the FDs of the low-z and unknown-z popu- 
lations are indistinguishable, having a KS probability of 0.42. 

That the SIDs and FDs both favor a relationship between the 
unknown-z and low-z HSP populations provides some confi- 
dence that this may, in fact, be the case. This is supported by 
the redshift distribution in the a-Fss plane, depicted in Figure 
10. HSPs with high redshifts are clustered at low fluxes and 
hard a. This is not unexpected given an upper limit upon the lu- 
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minosity of HSPs. The dotted lines in Figure 10 show constant 
redshift curves in the F^^-a plane for a lOOMeV-lOOGeV lu- 
minosity of 2 X 10"*^ ergs"' (the definition of in Abdo et al. 
2010c, note the distinction with the definition of F35). If all 
HSPs have intrinsic luminosities less than 2 x lO'^^'ergs"', no 
sources at a given z should be found to the right of the asso- 
ciated line. Since the volume of the visible Universe is dom- 
inated by z ~ 1, the majority of high-z objects should then be 
found up against the instrumental flux limit, i.e., at low fluxes 
and nearly flat spectra. Sources with harder spectra are likely 
to have higher bolometric gamma-ray luminosities since they 
are likely to be more below the inverse-Compton peak, biasing 
high-z sources towards harder spectra. 

However, the HSPs without redshifts are noticeably absent 
within this high-z dominated region. This is responsible for 
the fact that they more closely share their SID and FD with 
the low-z HSPs. The presence of HSPs without redshifts at 
a variety of as and F35S suggests that this is not a result of a 
selection effect upon either It is nonetheless possible that there 
is some instrumental effect which prevents the measurement of 
high redshifts in intrinsically bright, soft objects, in which case 
a population of high-luminosity HSPs with high (and therefore 
unmeasured) zs would necessarily be located in regions with 
predominantly low measured zs. However, this is belied by the 
broad distribution of luminosities and spectral indexes among 
the HSPs with measured redshifts. Thus, we adopt the simpler 
explanation: the HSPs with unknown redshifts have and SID 
and FD similar to low-z objects because they are intrinsically 
dim, nearby objects. 
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